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Abstract 



In this PhD thesis, we develop models of surfaces of epitaxially grown crystals. In the in- 
troductory chapter ||, we introduce the most important physical processes which occur on 
crystal surfaces in a molecular beam epitaxy (MBE) environment. In the first part of this 
work (chapters |2|-Q) we model the (001) surface of the II- VI semiconductors CdTe and ZnSe, 
our main focus being on CdTe. In the second part (chapters ^ and we study generic fea- 
tures of epitaxial growth which are not specific to one particular material. These are kinetic 
roughening and the formation of mounds. Both effects can be investigated only after the 
deposition of a thick film on the substrate. Therefore, we must restrict ourselves to simple 
models which can be simulated with moderate computational effort. 

In chapter ^, we introduce a lattice gas model of flat (001) surfaces of CdTe and ZnSe 
in thermal equilibrium. Both the arrangement of metal atoms in vacancy structures and 
the dimerization of nonmetal atoms are considered. The surface is represented by a two- 
dimensional square lattice. We introduce anisotropic interactions between nearest and next 
nearest neighbour sites which reflect the known properties of the surface. The phase diagram 
of this model is determined by means of transfer matrix calculations and Monte Carlo simu- 
lations. The phases of the model are identified with different reconstructions which have been 
observed in experiments. In particular, the transition between a Cd terminated c(2 x 2) recon- 
struction at low temperature and a Cd terminated (2 x 1) reconstruction at high temperature 
which occurs on the CdTe surface can be explained as an encompanying effect of an order- 
disorder phase transition. Here, the small energy difference between both reconstructions 
plays an important role. 

Crystal surfaces in an MBE environment are not in thermal equilibrium. Instead, non- 
equilibrium processes like growth and sublimation occur. Are the essential features of our 
equilibrium model preserved under these conditions or are there completely different, nonequi- 
librium effects? To answer this question, we need to extend our two-dimensional lattice gas 
to a model of a three-dimensional crystal which can be used to simulate growth and subli- 
mation. The basic idea is that there is an isotropic attraction between particles in the bulk 
of the crystal while particles on the surface interact with the anisotropic interactions of the 
lattice gas. 

Such a model which, however, contains several simplifications compared to a realistic 
model of CdTe, is investigated in chapter |3|. We neglect the dimerization of surface Te atoms 
and simulate a simple cubic lattice instead of the zinc-blende lattice of CdTe. In spite of these 
simplifications, many effects which have been observed on sublimating CdTe(OOl) surfaces in 
vacuum and under an external flux of pure Cd or pure Te can be reproduced. Additionally, the 
behaviour of this model agrees qualitatively with that of a simplified version of our lattice gas 
where the dimerization of Te atoms is neglected. Thus, we have shown that the investigation 
of two-dimensional lattice gas models in thermal equilibrium is an appropriate tool to obtain 
insight into the behaviour of reconstructed surfaces in an MBE environment. On the other 
hand, the non-equilibrium conditions of sublimation induce characteristic deviations from the 
behaviour of the equilibrium model which should be considered in precise investigations. 

Starting from this result, we develop a model of the CdTe(OOl) surface which is closer to 
reality. It is introduced in section ^. We simulate a zinc-blende lattice instead of a simple 
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cubic lattice. However, we still neglect the dimerization of Te atoms. An appropriate choice 
of the parameter set yields a rough quantitative agreement with experimental results. This 
indicates, that an extension towards a truly realistic model of CdTe(OOl) should be possible 
as soon as experimental data are available which allow for an adaptation of the model. We 
perform simulations of atomic layer epitaxy (ALE). Although this process is widely used in 
technological applications and basic research, to our knowledge it has never been investigated 
by means of Monte Carlo simulations before. We find a growth rate of one half monolayer 
per cycle and an alternation between a rough and a flat surface in subsequent cycles. Both 
effects have also been observed in experiments. 

In chapter ^, we investigate the fractal properties of kinetically roughened surfaces. We 
consider a model of an unreconstructed surface of a simple cubic crystal during MBE growth 
at low temperature. This model yields self-similar, fractal surfaces. We apply the WTMM 
method to obtain the first complete singularity spectra of such surfaces. The statistical 
properties of growing surfaces are invariant under dynamic scaling. We identify the dynamic 
exponent a with the Holder exponent which maximizes the singularity spectrum. This scaling 
behaviour is a generalization of the well-known Family- Vicsek scaling. However, we measure 
a strong dependency of the scaling exponents on the desorption rate which indicates that they 
are non-universal. 

Finally, in chapter ^ we consider the coarsening of mounded surfaces which are created 
during growth if the diffusion of particles across steps is suppressed by a strong Schwoebel 
barrier. A continuum model of this process which was studied first by Siegert et. al. predicts a 
dependency of the dynamic exponents on the symmetry of the surface. To test this prediction, 
we simulate growth in the presence of an infinite Schwoebel barrier on different lattice struc- 
tures. We use an efficient simulation algorithm which traces the motion of a single particle 
from deposition until it has reached a state where it is strongly bound. We find, contrary to 
the results of the continuum theory, that the scaling exponents are independent of the surface 
symmetry. Instead, the exponents depend on the speed of the diffusion of particles along step 
edges. The exponents which are measured in the limit of strong step edge diffusion can be 
explained by means of a simple consideration. 

In the two appendices we explain the techniques which we have used to investigate the 
models which are considered in this thesis. These are the transfer matrix method and contin- 
uous time Monte Carlo simulations. 
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Zusammenfassung 



In dieser Doktorarbeit werden Modelle fiir die Oberflachen epitaktisch gewachsener Kristalle 
entwickelt. Im Einleitungskapitel || werden die wichtigsten Prozesse auf Kristalloberflachen in 
einer fiir die Molekularstrahlepitaxie (MBE) typischen Umgebung vorgestellt. Im ersten Teil 
der Arbeit (Kapitel § prasentieren wir Modelle fiir die (OOl)-Oberflachen der II- VI Halb- 
leiter CdTe und ZnSe. Hauptsachlich beschaftigen wir uns dabei mit CdTe. Im zweiten Teil 
(Kapitel || und |6|) werden allgemeine Eigenschaften des Wachstumsprozesses untersucht, die 
nicht an ein spezielles Material gebunden sind. Dies sind das kinetische Aufrauen wachsender 
Oberflachen und das Hiigelwachstum. Da beide Effekte nur nach dem Aufdampfen dicker 
Schichten beobaclitet werden konnen, miissen wir uns bei ihrer Untersuchung auf einfache 
Modelle beschranken, die mit relativ wenig Rechenaufwand simuliert werden konnen. 

In Kapitel || stellen wir ein Gittergasmodell fiir flache CdTe(OOl)- und ZnSe(001)-Oberfla- 
chen im thermischen Gleichgewicht vor. Sowohl die Anordnung der Metallatome in Leerstel- 
lenstrukturen als auch die Dimerisierung der Nichtmetallatome werden beriicksichtigt. Die 
Oberflache wird durch ein zweidimensionales Quadratgitter dargestellt. Wir fiihren anisotrope 
Wechselwirkungen zwischen nachsten und iibernachsten Nachbarplatzen ein, in denen die be- 
kannten Eigenschaften der Oberflache beriicksichtigt werden. Das Phasendiagramm des Mo- 
dells wird mit Transfermatrixrechnungen und Monte Carlo Simulationen bestimmt. Die ein- 
zelnen Phasen werden mit verschiedenen experimentell untersuchten Rekonstruktionen iden- 
tifiziert. Insbesondere kann der Ubergang von einer Cd-terminierten c(2 x 2)-Rekonstruktion 
bei tiefen Temperaturen zu einer ebenfalls Cd-terminierten (2 x l)-Rekonstruktion bei ho- 
hen Temperaturen, der auf CdTe-Oberflachen beobachtet wurde, als Folge eines Ordnungs- 
Unordnungs-Phaseniibergangs erklart werden. Dabei spielt der kleine Energieunterschied zwi- 
schen den beiden Rekonstruktionen eine wesentliche Rolle. 

Kristalloberflachen in einer MBE-Anlage befinden sich nicht im thermischen Gleichge- 
wicht. Stattdessen finden Nicht gleichgewichtsprozesse wie Wachstum und Sublimation statt. 
Bleiben die wesentlichen Eigenschaften unseres Gleichgewichtsmodells unter diesen Bedin- 
gungen erhalten oder treten vollkommen neue Effekte auf? Um diese Frage zu beantworten, 
miissen wir unser zweidimensionales Gittergas zu einem Modell eines dreidimensionalen Kris- 
tails erweitern, mit dem Wachstum und Sublimation simuliert werden konnen. Die grund- 
legende Idee dabei ist, dass sich Atome im Inneren des Kristalls isotrop anziehen, wahrend 
zwischen Teilchen an der Oberflache die anisotropen Wechselwirkungen des Gittergases wir- 
ken. 

Ein solches Modell, das im Vergleich zu einer realistischen Simulation des CdTe allerdings 
stark vereinfacht ist, betrachten wir in Kapitel ^. Wir vernachlassigen die Dimerisierung 
von Te-Atomen an der Oberflache und simulieren ein einfach kubisches Gitter anstelle des 
Zinkblendengitters des CdTe. Trotz dieser Vereinfachungen konnen mit diesem Modell viele 
Effekte reproduziert werden, die auf sublimierenden CdTe(001)-Oberflachen im Vakuum und 
unter einem externen Fluss aus reinem Cd oder reinem Te beobachtet wurden. Weiterhin 
stimmt das Verhalten des Modells qualitativ mit dem einer vereinfachten Version unseres Git- 
tergases iiberein, in dem die Dimerisierung von Te-Atomen nicht beriicksicht wird. Damit ist 
gezeigt, dass die Untersuchung zweidimensionaler Gittergasmodelle im thermischen Gleichge- 
wicht prinzipiell ein geeignetes Werkzeug zum Verstandnis rekonstruierter Oberflachen unter 
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MBE-Bedingungen ist. Andererseits verursachen die Nichgleichgewichtsbedingungen der Sub- 
limation aber auch charakteristische Abweichungen vom Verhalten des Gleichgewichtsmodells, 
die in genauen Untersuchungen beriicksichtigt werden sollten. 

Ausgehend von diesem Ergebnis entwickeln wir ein realitatsnaheres Modell der CdTe(OOl)- 
Oberflache, das in Kapitel |^ vorgestellt wird. Wir simulieren das Zinkblendengitter anstelle 
eines einfach kubischen Gitters. Die Dimerisierung von Te-Atomen wird aber weiterhin ver- 
nachlassigt. Durch Wahl eines geeigneten Parametersatzes konnen experimentelle Ergebnisse 
nicht nur rein qualitativ, sondern auch grob quantitativ reproduziert werden. Das deutet 
darauf hin, dass eine Erweiterung zu einem wirklich realistischen Modell moglich sein sollte, 
sobald geeignete experimentelle Daten zur Anpassung vorliegen. Mit diesem Modell werden 
Simulationen der Atomlagenepitaxie (ALE) durchgefiihrt. Obwohl dieser Prozess in Tech- 
nik und Grundlagenforschung haufig angewendet wird, wurde er unseres Wissens nach noch 
nie in Monte Carlo Simulationen untersucht. Wir beobacliten eine Wachstumsrate von einer 
halben Monolage pro Zyklus und einen Wechsel zwischen rauen und glatten Oberflachen in 
aufeinanderfolgenden Zyklen. Beide EfFekte wurden auch in Experimenten beobachtet. 

In Kapitel |5| werden die fraktalen Eigenschaften kinetisch aufgerauter Oberflachen un- 
tersucht. Wir betrachten ein Modell einer unrekonstruierten Oberflache eines einfach ku- 
bischen Kristalls wahrend des MBE-Wachstums bei niedriger Temperatur. Dabei entste- 
hen selbstahnliche Oberflachen, deren multifraktale Spektren mit Hilfe der WTMM-Methode 
erstmals vollstandig bestimmt werden konnten. Die statistischen Eigenschaften wachsender 
Oberflachen sind invariant unter der dynamischen Skalentransformation. Wir identifizieren 
den dynamischen Exponenten a. mit demjenigen Holder-Exponenten, der das multifraktale 
Spektrum maximiert. Dieses Skalenverhalten stellt eine Verallgemeinerung der aus der Li- 
teratur bekannten Family- Vicsek-Skaleninvarianz dar. Allerdings messen wir eine starke 
Abhangigkeit der Skalenexponenten von der Desorptionsrate, was darauf hindeutet dass diese 
nicht universell sind. 

Abschliefiend untersuchen wir in Kapitel ^ die Vergroberungsdynamik hiigeliger Ober- 
flachen, wie sie wahrend des Wachstums entstehen, wenn die Teilchendiffusion iiber Stufen 
hinweg durch eine starke Schwobelbarriere unterdriickt wird. Ein Kontinuumsmodell dieses 
Vorgangs, das zuerst von Siegert et. al. untersucht wurde, sagt eine Abhangigkeit der dyna- 
mischen Skalenexponenten von der Symmetric der Oberflache voraus. Um diese Vorhersage 
zu liberpriifen, simulieren wir das Wachstum mit unendlich hoher Schwobelbarriere auf ver- 
schiedenen Gitterstrukturen. Dabei wird ein effizienter Simulationsalgorithmus verwendet, in 
dem jeweils die Bewegungen eines einzelnen Teilchens von der Deposition bis zum Erreichen 
eines fest gebundenen Endzustands verfolgt werden. Wir stellen fest, dass die Skalenexponen- 
ten, anders als von der Kontinuumstheorie vorhergesagt, nicht von der Oberflachensymmetrie 
abhangen. Stattdessen beobachten wir eine Abhangigkeit der Exponenten von der Geschwin- 
digkeit der Diffusion entlang von Stufenkanten. Die bei starker Kantendiffusion gemessenen 
Skalenexponenten konnen mit einer einfachen Abschatzung erklart werden. 

In den beiden Anhangen werden die Techniken, mit denen die in dieser Arbeit betrachteten 
Modelle untersucht wurden, also die Transfermatrixmethode und die Monte Carlo Simulation 
mit kontinuierlicher Zeit, vorgestellt. 
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Chapter 1 



The physics of crystal surfaces in an 
MBE environment 



In the fabrication of modern semiconductor devices like computer chips or laser diodes one 
frequently is interested in depositing thin films on a monocrystalline substrate. Important 
techniques to achieve this goal are molecular beam epitaxy (MBE) and related methods |44]. 



A sketch of an MBE chamber is shown in figure LI. The whole process is performed in ultra 
high vacuum (UHV) to avoid the intrinsic instabilities of crystal growth in a liquid or gassy 
medium which result from diffusive and convective processes . Material is evaporated in an 
effusion cell and creates an atomic or molecular beam which is directed towards the substrate. 
The mean free path of the particles in the beam is large compared to the distance between the 
effusion cell and the substrate such that the particles move on straight lines. Thus, the whole 
surface is exposed to an even material flux. The success of MBE depends on the ability to 
control the physical processes which occur on the surface of the substrate itself. This might 
be done e.g. by a regulation of temperature and particle flux. 

UHV chamber 




Figure 1.1: Schematical sketch of an MBE 
growth system. In an UHV chamber, material 
which is evaporated in an effusion cell creates 
an atomic or molecular beam. Particles from 
the beam are attached to the substrate and cre- 
ate a film. 



MBE growth leads to the formation of a rich variety of surface morphologies - to the regret 
of many practicians who are interested in smooth films. From the theoretical point of view, 
these phenomena pose a challenge to our understanding. First, the physics of crystal surfaces 
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is determined by a complex interplay of various effects like surface reconstruction and the 
adsorption, diffusion and desorption of adatoms. Second, a crystal under MBE conditions is 
in a state far from thermal equilibrium. In this case, no general theoretical framework like 
equilibrium statistical mechanics is available. Instead, the behaviour of the system might be 
determined by the details of microscopic processes which occur on atomic lengthscales and 
short timescales of 10~^^ s. However, large surface formations like mounds have extensions of 
up to lfj,m and their formation takes several minutes to hours. This wide range of relevant 
length- and timescales has inspired a variety of approaches which range from first principles 
methods based on quantum mechanics to macroscopic differential equations. In this thesis, 
we study lattice gas models which build a bridge between a purely microscopic and a purely 
macroscopic point of view. In contrast to differential equations, they consider the atomistic 
nature of matter. On the other hand, their conceptual simplicity permits the investigation of 
large systems on long timescales, which is not possible by means of molecular dynamics or 
density functional methods. A large fraction of our work is devoted to II- VI semiconductors, 
in particular CdTe and ZnSe. However, we also investigate more abstract models which aim 
at an understanding of general, material-independent features. 

In the description of a non-equilibrium process like growth or sublimation of a crystal, 
two different aspects have to be considered. On the one hand, systems which are not in 
thermal equilibrium tend to approach to it. Therefore, an understanding of the evolution of 
such a system requires an understanding of its equilibrium configuration. In this thesis, we 
present the results of models of pure surfaces, i.e. surfaces without adsorbates. Here, we have 
to consider surface reconstructions and the roughening transition. On the other hand, the 
way the system evolves depends on the kinetics of processes which occur on the surface like 
deposition, diffusion and desorption of atoms and molecules. 



1.1 Surface reconstruction 

In the groundstate, i.e. at temperature T = 0, atoms in a crystal arrange in such a way that 
the energy of the system is minimized. This principle determines the material-specific lattice 
structure. Atoms at the surface of the crystal have a smaller number of binding partners than 
atoms in the bulk. Therefore, it is frequently energetically favourable if the surface atoms 
arrange in a structure which has a lower degree of symmetry than a planar section of the 
bulk. This phenomenon is called surface reconstruction. 

Reconstructions have been observed in various materials. In II- VI semiconductors, a 
variety of reconstructions consisting of vacancy structures, dimers and trimers has been found 
pl| . These will be discussed in detail in chapers ^-Q where we present models of reconstructed 
(001) surfaces of CdTe and ZnSe. 

In general, a theoretical prediction whether a surface of a particular material reconstructs 
and which reconstructions can be found is a hard problem even at T = 0. It requires an 
investigation of the properties of the chemical bonds in the crystal. For this purpose, molecular 



statics using empirical potentials [99| and density functional theory (DFT) |37, ^ have 
been employed. The computational burden of these methods is comparatively high, which 
restricts their practical applicability to systems consisting of only a few atoms. Due to the 
periodicity of crystal surfaces this is not a severe restriction if one is interested in ground state 
properties. 

In the case of polar semiconductors like III-V and II- VI compounds, there is a simple 
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consideration which allows to exclude energetically unfavourable surface configurations fl^ , 
8C]. Consider the electronic structure of a compound semiconductor. A schematical sketch of 
the energy levels in ZnSe is shown in figure |L^. Other compound semiconductors like GaAs 
and CdTe show a similar behaviour. Prom the s- and p-orbitals of the valence electrons sp^ 
hybrids are created which have an energy e^yb . The energy bands of the crystal are obtained 
from linear combinations of these atomic orbitals. The bonding combinations yield the valence 
band, the antibonding ones the conduction band. At T = 0, the valence band is fully occupied 
and the conduction band is empty. 
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Figure 1.2: Electronic states in ZnSe 
according to [^. All energy levels are 
given in eV. Eg and £p are the ener- 
gies of the s- and p-orbitals, respec- 
tively. These hybridize to form sp^ or- 
bitals which have an energy Shyb- The 
dangling bond energy of anions is inside 
the valence band, that of cations is in 
the conduction band. 



In the bulk of the crystal, the four sp^ orbitals of each atom overlap with orbitals of 
neighbours and create chemical bonds. Each bond is occupied by two electrons. Due to the 
different electronegativities of the elements the metal atoms carry positive charge (cations) 
and the atoms of the nonmetallic elements are negatively charged (anions). Atoms at the 
surface have a smaller number of nearest neighbours such that some of their orbitals are 
directed out of the crystal. These are denoted as dangling bonds. Cation dangling bonds are 
empty since their energy level is in the conduction band. On the contrary, the binding energy 
of electrons in the anion dangling bonds is below the valence band maximum such that each 
of these orbitals is occupied by two electrons. The total charge of each atom is the sum of its 
valence nuclear charge, the charge of the bonding electrons and the charge of the electrons in 
the dangling bonds. In general, it is nonzero. 

In the majority of all conceivable arrangements of atoms in the terminating layer of the 
crystal there is a net charge of the surface which creates an electric field in the bulk. Energeti- 
cally, this is extremely unfavourable [^] . Consequently, in nature only surface reconstructions 
without a net charge are realized. This condition is fulfilled if the number of electrons is just 
sufficient to fill the chemical bonds and the anion dangling bonds. In the literature, this law 



is known as electron counting rule |8C]. 

At higher temperature, the properties of the surface are influenced by thermodynamic 
effects. Frequently, phase transitions between different reconstructions are observed. Strictly 
speaking, phase transitions exist only in the thermodynamic limit of an infinite number of 
degrees of freedom. Therefore, their theoretical investigation requires the study of systems 
with a large number of atoms which is beyond the scope of first principles methods or molecular 
dynamics simulations using realistic empirical potentials. In this case, effective models which 
neglect several details but preserve the relevant features of the atomic interactions are required. 
The models of CdTe and ZnSe which we present in chapters belong to this class. 
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1.2 The roughening transition 



In thermal equilibrium at low temperature, surfaces of homogeneous crystals which are ori- 
ented in a high symmetry direction are essentially flat apart from the pattern of a possible 
reconstruction. Since atoms on the surface are bound more weakly than atoms in the bulk, the 
system minimizes its energy by reducing the area of the surface. The most frequent thermal 
excitations are single atoms or small groups of atoms diffusing around on the surface and the 
vacancies they have left when jumping out of the terrace. 

The creation of a step of atomic height and length / costs an energy estcp^- In good 
approximation, the step energy per unit length Egtep is independent of temperature. Since 
atoms at the step edge have a smaller number of binding partners than atoms which are 
incorporated in a terrace, estcp is positive. However, at the creation of the step an entropy 
Sstcp^ is gained since there are many step configurations with equal energy. Consequently, the 
free energy per unit length, 7 = Egtep — Tsgtcp decreases with increasing temperature. 7 is 
denoted as line tension of steps. At a critical temperature Tr, there is a Kosterlitz-Thouless 
phase transition |^, At higher temperature, we have 7 = 0. Then, steps are created 
spontaneously by thermal fluctuations and the surface becomes rough. This phase transition 
is called roughening transition. 

A quantitative measure of the roughness of a surface is the height-height correlation func- 
tion G{1). Consider a system of cartesian coordinates {x,y,z) the z-axis of which is per- 
pendicular to the surface. In the remainder of this thesis we will use this orientation of the 
coordinate system. We define the height of the surface h{x) as function of the vectors x in the 
two-dimensional {x, y) space which returns the z-coordinate of the surface at x. Then, G{1) 
is defined as 

G{l):=(^(^hix)-h{x + l)y^ . (1.1) 

Below Tr, one finds that G{1) remains finite in the limit of an infinite distance |/|. This 
expresses the fact that the surface is essentially flat. On the contrary, at T > Tr we have 
G{1) ~ ln(|/|). This logarithmic divergence corresponds to a surface with strong fluctuations. 

Usually, molecular beam epitaxy is performed at temperatures far below the roughening 
transition. This is also the case in the growth models which are investigated in this thesis. 
Consequently, surfaces of homogeneous crystals in thermal equilibrium should be flat at typical 
MBE temperatures. 



1.3 Strained overlayers 

In heteroepitaxy, the deposited fllm consists of a different material than the substrate. If 
substrate and adsorbate have different lattice constants and a^, equilibrium conflgurations 
where the surface is not flat can be found even at low temperature. At the adsorbate-substrate 
interface, the interaction of the adsorbate particles with the substrate forces the adatoms to 
adapt to the substrate lattice constant a^- However, the adatom-adatom interaction favours a 
distance between adsorbate particles. This induces strain in the adsorbate. The equilibrium 
configuration depends on the relative strength of the competing forces and the relative lattice 
mismatch (cq — as)/as. 

Three different growth modes have been observed |82, If the attraction between 

adsorbate particles is much stronger than the adsorbate-substrate interaction and/or the 
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lattice mismatch is large, one obtains small droplets of the adsorbate. This is called the 
Volmer- Weber growth mode. Conversely, if the binding of the adsorbate atoms to the substrate 
is stronger than the adsorbate-adsorbate bond, small amounts of the adsorbate form a flat 
layer. The behaviour of a thicker adsorbate layer depends on whether the adsorbate-adsorbate 
interaction is strong enough to outweigh the elastic energy of the strain. If this is the case, 
thick layers are flat. This is the Frank-van der Merwe mode. Otherwise, in the Stranski- 
Krastanov mode, adsorbate islands are formed on top of a thin wetting layer. In very thick 
films, dislocations are created. Their creation costs energy, however, they permit the adsorbate 
to adapt to its own lattice constant a^. Recently, specialized simulation techniques have been 
developed which allow for the simulation of strained heteroepitaxial growth |70, |7l|, They 
are computationally expensive, which restricts their practical applicability to 1+1-dimensional 
models. 

In the remainder of this thesis we restrict ourselves to homoepitaxial growth where sub- 
strate and adsorbate particles are identical. 



1.4 Motion of atoms on crystal surfaces 

Simulations of surfaces in an MBE environment are based on rules which describe the motion 
of atoms. In growth models, the atomic nuclei are treated as classical particles the motion of 
which is, in principle, described by Newton's laws. Of course, the forces which act between 
the particles are determined by the electronic structure which follows quantum-mechanical 
laws. In molecular dynamics (MD) simulations, the equation of motion of the particles is 
solved numerically. The forces are calculated by means of empirical potentials or DFT. 

In the energy landscape of a large number of atoms in phase space, there are many local 
minima. These correspond to configurations where the atoms arrange e.g. in a regular crystal 
lattice. If the heights of the barriers between the minima are noticeably higher than kT, most 
of the time the actual state of the system will fluctuate around one of them. These are the 
thermally excited lattice vibrations of the crystal. Very rarely, the system leaves a minimum 
and jumps into a neighbouring one. These are processes where the arrangement of atoms 
changes like e.g. a diffusion jump of an adatom on the surface. In a MD simulation, the lattice 
vibrations are considered explicitly. This is the reason why this method is computationally 
extremely expensive. On modern computers, the motion of several thousand atoms can be 
traced for at most some nanoseconds, which is by far too short to investigate MBE. 

In Monte Carlo simulations, we investigate the system on a coarse grained scale. Only 
changes of the arrangement of atoms are considered explicitly, while the chaotic motion of 
atoms due to lattice vibrations is considered implicitly in a statistical manner. It determines 
the rate p of events where the atomic arrangement changes. One can show quite generally 
1 ^ that for thermally excited events p is given in good approximation by an Arrhenius law 

E is the height of the energy barrier which must be overcome to change from the initial to 
the final minimum, z/ is an attempt frequency which is on the order of magnitude of the 
Debye frequency of the crystal. In general, it is reasonable to assume that it is identical for all 
processes which occur on the surface. Compared to MD, the Monte Carlo technique reduces 
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the computational effort by many orders of magnitude and makes it feasible to consider the 
typical timescales of MBE. 

Frequently, the position of the atoms in the relevant local energy minima can be assigned to 
sites in a (potentially distorted) crystal lattice. In this case, lattice gas models are appropriate. 
We consider a fixed crystal lattice. Each lattice site is either occupied by an atom or empty. 
At the deposition of an atom, the state of one site changes from empty to occupied. If an atom 
is desorbed, its site is switched from occupied to empty. At the diffusion of an atom, the states 
of two sites are exchanged. In the presence of surface reconstructions, the topology of the 
lattice at the surface might be different from that in the bulk. It is possible to extend a lattice 
gas model to include this effect by adding additional discrete degrees of freedom. Examples 



are the formation of Si dimers on Si(OOl) in the simulations of Rockett et. al. |^, 94| and the 
dimerization of Te atoms in our model of the CdTe(OOl) surface which will be presented in 
section ^. 

In addition to the application in Monte Carlo simulations, lattice gases and Arrhenius 
activated processes are the basis of most theories which consider the dynamics of processes 
which occur at crystal surfaces. Examples include nucleation theory |122| , |128| and BCF 

1.4.1 Deposition 

Atoms impinge on the substrate with a kinetic energy ~ kT^s , where Tes is the temperature of 
the effusion cell. Typically, Tcfj is higher than the substrate temperature. Additional energy 
is gained from the formation of chemical bonds to the surface. Consequently, the energy of the 
newly arriving particles is considerably higher than that of particles on the surface which are 
already thermalized at substrate temperature. This leads to an increased mobility of adatoms 
until the excess energy is dissipated. The particles get trapped at energetically favourable sites 
with high coordination number. Clearly, such sites will be found preferentially in layers close 
to the substrate. Therefore, the motion of the adatoms is biased in the "downhill" direction. 
The momentum of the particles in the beam which is directed towards the substrate enhances 



this bias. In the literature |38, 13C] results of molecular dynamics simulations have been 
reported which confirm these qualitative considerations. 

In our simulations, we consider this effect by means of effective rules termed incorporation 
and downhill funneling. In incorporation, a particle moves to the lattice site at the lowest 
height within an incorporation radius ri around the position where it hits the surface |103]. 



The simulations presented in [l38| , |130|| suggest that reasonable values of rj should be on the 
order of magnitude of a few lattice constants. In downhill funneling |32, 33], adatoms perform 



a biased random walk to nearest neighbour sites lying lower than the current site. The particle 
stops when the height of all neighbour sites is equal to or greater than the current height. 

At the deposition of particles, precursor states can play an important role. Particles which 
land on the surface first enter a state where they are only weakly bound to the surface. Later, 
a strong chemical bond is formed. Meanwhile, they may diffuse over a considerable distance. 
In chapter we investigate a model where such weakly bound states play an important role. 

1.4.2 Difrusion 

To understand the properties of diffusion, the concept of the potential energy surface (PES) 
or surface energy map is employed. It is assumed that the energy barriers at zero temperature 
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are a reasonable approximation for the energy barriers at growth temperature. Consider an 
atom close to a surface. We fix its position in the x- and y-direction at x = {x, y). There is no 
constraint in the z-coordinate of this test atom and the positions of the atoms in the crystal. 
We determine the equilibrium configuration at T = by minimizing the total energy with 
respect to these remaining degrees of freedom. Measuring the total energy of the system as a 
function of x, we obtain the PES of an adatom. Local minima of the PES correspond to sites 
where adatoms are stable with respect to small perturbations e.g. from lattice vibrations. 

In the investigation of the hopping of adatoms by means of Monte Carlo simulations, we 
consider processes where atoms move between these sites. The energy barrier that a single 
atom must overcome to hop to a neighbouring site equals the difference between the energy 
of the initial configuration and the saddle point of the PES which lies between the initial 
and the final site. This information is sufficient to calculate the ratios between the rates of 



various diffusion processes which occur on the surface from the Arrhenius law (equation 1.2). 
Interactions between atoms are short-ranged. Consequently, the rates of diffusion processes 
depend only on the local environment of the diffusing atom within a range of a few lattice 
constants such that it is sufficient to consider a limited number of typical configurations. 




Figure 1.3: Potential energy 
of a test atom on the surface 
of a two-dimensional Lennard- 
Jones crystal at zero tempera- 
ture function of its posi- 
tion. By courtesy of F. Much. 



At typical MBE growth temperatures, bulk diffusion can be neglected. The dominant 
mechanism of material transport is surface diffusion. We discuss some important effects in 
the simplified case of the surface of a two-dimensional crystal of Lennard-Jones particles 
||70|, [7l|. This is an example of a system without surface reconstruction. Its PES is shown in 



figure 1.3. As expected, the potential energy of the test particle has local minima at lattice 
sites. Therefore, at the deposition of atoms the crystal grows in its own lattice structure such 
that a lattice gas modelling of this system would be appropriate. At a step edge, where it is 
surrounded by a large number of neighbours an atom is bound stronger than on a flat surface. 
Consequently, there is a high energy barrier for the detachment of an atom from the step 
edge. Similarly, adatoms are bound to other adatoms. 

An adatom wich crosses a step edge from above must overcome a greater energy barrier 
than that for diffusion on a flat terrace. This so-called Schwoebel effect |106, 107] originates 
from the fact that a particle at the rim of the step has only a small number of binding partners. 
Therefore, hopping down a step is suppressed significantly compared to terrace diffusion. 

The same effects are found in the physically relevant case of the surface of a three- 
dimensional crystal. Additionally, we have to consider the diffusion of atoms along step 
edges. In many materials, step edge diffusion at straight steps is as fast as diffusion on a flat 
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terrace or even faster. Atoms in kinks of an edge have a high coordination number. Con- 
sequently, particles get trapped at these kink sites. Similar to the Schwoebel effect, energy 
barriers which suppress step edge diffusion around corners are possible. 

In addition to the hopping of atoms, diffusion can take place via exchange processes. An 
adatom pushes another atom out of the surface and occupies its site. The atom pushed out 
becomes an adatom at a neighbouring site. From the viewpoint of Monte Carlo simulations 
of homoepitaxy both processes are equivalent, since initial and final state of hopping and 
exchange diffusion are identical. However, the energy barriers and therefore the rates of both 
processes are different. Whether diffusion occurs preferentially via hopping or exchange is 
material dependent. 

The PES of a reconstructed surface has the symmetry of the reconstruction. This may 
lead to strong anisotropics in the diffusion rates. 

In the literature, PES and diffusion barriers for various materials have been published 



which have been calculated by means of empirical potentials [^, 97, 127 1 and density 



functional theory pq, 54, 55, 115 1. However, currently no such data are available for the 



materials which we investigate in this thesis. Therefore, in our models we will use param- 
eterizations of energy barriers which consider the effects which have been discussed in this 
section. 



1.4.3 Desorption 

Similar to diffusion, the desorption of atoms from the surface is a thermally activated process 
the rate of which is given by an Arrhenius law. In general, the higher the coordination number 
of an atom, the greater is the energy barrier for desorption. Consequently, the desorption of 
single adatoms which diffuse on a flat terrace occurs at the highest rate, while atoms which 
are attached to step edges or incorporated in a flat surface are much more stable. In most 
cases it is reasonable to assume that the energy barrier which must be overcome to remove 
an atom from the surface equals the difference between the energy of the surface without the 
adatom and a single atom in vacuum and the energy of the crystal with the adatom at the 
surface. 



1.5 Epitaxial crystal growth 

In the following, we measure lengthscales in appropriate units such that the rate at which 
particles are deposited at one site equals the particle flux F, while diffusion constants of 
adatoms are identical to hopping rates. This can be done by setting the unit length to a 
small multiple of the lattice constant of the substrate which depends on the geometry of the 
crystal lattice. The average time an atom stays on the surface before it is incorporated into 
the crystal is 1/F. Diffusion and desorption processes which occur on much longer timescales 
can be neglected. Typically, this is the case for the detachment of particles from step edges 
and for the decay of islands which consist of more than i* atoms. The critical nucleus size i* 
is small at typical MBE temperatures. Frequently, even pairs of adatoms are stable such that 
i* = 1. Adatoms which are deposited on the surface diffuse around until they stick to a step 
edge or collide with more than i* other adatoms and form a stable island. The latter process 
is called nucleation. Once an island is formed, it captures more and more adatoms and grows. 
Thereby, the adatom density in the vicinity of the island is reduced which suppresses further 
nucleation events. Thus, islands keep a typical distance li to their neighbours. 
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Nucleation theory [122, 128] states that the minimal value of li which is obtained during 
the deposition of the first monolayer on a perfectly flat surface has a powerlaw dependence 
on the ratio of the adatom diffusion constant D and F: 

li oc {^^ (1.3) 

The value of k depends on various parameters like i* , the fractal dimension of islands, the 
mobility of small clusters of adatoms and on whether a significant fraction of adatoms is 
desorbed. In the simplest case of i* = 1, compact islands and vanishing desorption, we have 



K = 1/6. We refer the reader to refs. |122, |128| | for a discussion of this subject. In the late 
phase of this suhmonolayer regime, neighbouring islands coalesce and li increases again. 

A perfectly flat surface is never realized in experiments although it is the equilibrium 
configuration at growth temperature. On the one hand, this is due to the fact that preparation 
techniques yield an intrinsic surface morphology and thermal smoothing is extremely slow. 
On the other hand, steps are created during growth by kinetic roughening and the Schwoebel 



instability which will be discussed in sections 1.5.1 and [1.5. 2. As long as the average distance 



Is between steps is much greater than li, the influence of the steps is weak and growth proceeds 
essentially the same way as on a flat surface. However, if Ig <^ h all adatoms are captured 
by the steps and no islands nucleate. In this case, growth proceeds in step flow mode. The 
influence of steps on growth can be investigated particularly well on a vicinal surface which 
consists of flat terraces separated by regularly spaced straight steps which are parallel to each 
other. Experimentally, vicinal surfaces can be created by cleaving a crystal in a small miscut 
angle to a high symmetry surface. 

1.5.1 Kinetic roughening 

Clearly, the deposition of material on a substrate in an MBE chamber is a stochastic process. 
Adatoms impinge on random positions. There are random time intervals between successive 
deposition events at one site. This randomness tends to create rough surfaces with a strongly 



fluctuating height function h{x) [|112| , 113]. On the other hand, diffusion tends to smooth the 
surface. In the following, we assume that growth starts at time t = 0. Since particles diffuse 
at a finite speed, diffusion can level out fluctuations only at lengthscales smaller than the 
correlation length ^(t) which increases with time. Properties of regions which are separated 
by a greater distance are uncorrelated. In particular, there is a random height difference 
between points at great distance. Consequently, growing surfaces are always rough at large 
lengthscales. This phenomenon is denoted as kinetic roughening, since it is induced by the 
kinetics of diffusing adatoms. In contrast to the roughening of surfaces in thermal equilibrium 
which occurs at high temperature, kinetic roughening is most effective at low temperature 
where diffusion is slow. 

The appearance of growing surfaces on small lengthscales depends on the properties of the 
diffusion process. In the absence of a Schwoebel barrier, surfaces with fractal properties can 
be found. In the modelling of these surfaces, frequently the concept of statistically self-affine 
surfaces is employed [^]. Mathematically, a self-affine surface is defined as a single- valued 
height function h{x) the properties of which are invariant under a simultaneous rescaling 
X — > 6x and h hb{h) ]|105|] . Here, b is an arbitrary real number. Clearly, an iterated 
transformation of the rescaled system with a scaling factor b' should yield the same result 
as a rescaling of the originial system with a factor bb'. This group property implies that 
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hy{hi){K)) = h}jy[h), which is fulfihed for hi){K) 
transformation is 

X bx; h 



b^h. Thus, the general form of the scale 

(1.4) 



b^h. 



The exponent H is denoted as Hurst exponent. In the special case H = 1, the surface 
is invariant under an isotropic rescaling of lateral extension and height. This property is 
denoted as self- similarity. 

In statistically self-affine surfaces, invariance under rescaling is not fulfilled for h{x) itself. 
Instead, we assume that there is a stochastic process which creates surfaces such that scale 
invariance is fulfilled for measurends averaged over an ensemble of independent realizations. 

Scale invariance implies conditions for the functional form of various quantities which 
can be used to measure the Hurst exponent. If x is scaled by a factor 6, the height-height 
correlation G{1) which was defined in equation |l.l| behaves like the square of the surface height 
such that G{bl) = b'^^G{l). This is fulfilled if and only if 



G{1) 




(1.5) 



where we have introduced I = Additionally, we consider the standard deviation of the 
surface height in a square of size L x L, 



w{L) 



dy{h{x)-{h{x))f 



(1.6) 



Under the scale transformation 1.4, a square of size L is transformed into a square of size bL. 
w{L) scales like the height function such that w{bL) = b^w{L), which implies 



w{L) = wqL". 



(1.7) 



Physical crystal surfaces are never statistically self-affine in the strict mathematical sense. 
First, there cannot be any structures which are smaller than the lattice constant a. This 



lower cutoff implies that invariance under the transformation L4 breaks down for too large 
b. Additionally, there is an upper cutoff due to the finite correlation length £,{t). At greater 
lengthscales, surface heights are uncorrelated which implies that w{L) = Woo = const, for 
L ^ ^(t) and G{1) = 2w%^ for I ^ ^{t). However, at late time where ^(t) is large, in several 



growth models there is a wide range of lengthscales a <C l,L <C ^(t) where equations |l.5| and 
|1.7| are fulfilled in good approximation. Surfaces with these properties are denoted as physical 
fractals. However, to fulfil symmetry properties like equation 1.4, the mean surface height 
must be subtracted. Therefore, we consider the reduced surface height f{x) := h{x) — {h{x)) 
instead of h{x) itself in the investigation of symmetry properties of growing surfaces. 

Careful investigations have shown, that some growth models yield surfaces with more 
complicated fractal properties which are denoted as multifractal |11, |2 
^ we will investigate such a model. 

1.5.2 Unstable growth 

On materials with a Schwoebel barrier, epitaxially grown surfaces have a mounded morphol- 
ogy. In the submonolayer regime, islands are formed at a typical distance li. Adatoms which 



24, 57]. In chapter 
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land on top of one of these islands diffuse around. Due to the Schwoebel barrier there is only 
a small probability that they cross the island edge and are incorporated in the layer below. 
Instead, they stay on top of the island until they either desorb or meet other adatoms and 
nucleate. This process creates islands in the second layer before the first layer is completed. 
Then, islands in the third layer nucleate on top of these islands and so on. Thus, in the course 
of time mounds build up. Their initial perimeter equals the distance between islands in the 
submonolayer regime. 

At the flank of a mound, there is a train of steps similar to a vicinal surface. Adatoms 
which approach a step edge from below are captured by the step. The attachment of atoms to 
steps from above is hindered by the Schwoebel effect. Nucleation on terraces can be neglected 
since the terrace width is necessarily smaller than li. Thus, the majority of the adatoms 
on a terrace is attached to the upper step. This yields a Schwoebel current in the uphill 
direction which tends to further increase the slope. However, processes like incorporation and 
downhill funneling keep the mounds from steepening infinitely. In these nonthermal processes, 
adatoms which impinge close to a step edge move in the downhill direction. This downhill 
current compensates the Schwoebel current at a magic slope rriQ at which the steepening of 
mounds stops. Frequently, neighbouring mounds coalesce such that the average perimeter of 
the mounds increases with time. This process is denoted as coarsening. 

We consider the initial phase of mound formation where the slope m is much smaller than 
mo in the case li ^ 1. Then, there are large terraces such that the influence of incorporation 
or downhill funneling can be neglected. In a simplifying manner we assume that the Schwoebel 
barrier is infinite such that there is no interlayer diffusion. In the first layer. Fa particles per 
unit time and unit area are deposited, where a is the fraction of the substrate which is not 
covered by the adsorbate. These particles cover the substrate such that 

^ = -Fa =^ a{t) = exp{-Ft). (1.8) 
dt 

Consequently, the first layer is not completed as long as m ^ niQ. This is the so-called Zeno 
effect [Q. The mounds are separated by deep grooves and there is no coarsening. In case of a 
finite Schwoebel barrier, there is an initial coarsening regime. However, coarsening stops at a 
critical mound size and deep grooves between the mounds are formed. Such a behaviour was 



predicted from continuous differential equations ^] and Monte Carlo simulations |25|, |4q]. 
Experimentally, it has been observed in homoepitaxial growth on Pt(lll) |4^ ]. 

When m approaches mg, incorporation processes become important. Then, the grooves are 
filled up and the flanks of the mounds assume their stable slope. Simultaneously, coarsening 
starts. In chapter 0, we will investigate this process. 



1.5.3 Dynamic scaling 

In the investigation of crystal growth, one frequently finds that surface snapshots are similar 
to rescaled surface snapshots taken at different time. In the growth of fractal surfaces, the 
correlation length ^(t) increases with time while the fractal properties do not change (figure 
[5.11 ). In the coarsening process of mounded surfaces, the size of the mounds grows but their 
shape remains the same (figure |6T| ). Therefore, it is plausible to assume that the statistical 
properties of the surface are invariant under a simultaneous transformation of spatial extension 
X, time t and the reduced surface height f{x,t) [|T^] : 

x^bx; t^bH; f^b"f. (1.9) 



19 



The powerlaw dependence of the scaHng factors for t and / on b is due to the group property of 
the transformation. This symmetry property is denoted as dynamic scaling. It is independent 



of the scale invariance 1.4 of self-affine surfaces. In many models of growth on an initially flat 
surface, it is valid after a transient phase like e.g. the initial formation of mounds in unstable 
growth. It is widely believed that the dynamic exponents a and z are universal. That means, 
that a and z are independent of "details" of the model, in particular the numerical values of 
parameters. Then, the values of a and z determine the universality class of a model. 

Dynamic scaling implies conditions for the time dependence of the surface width W{t) := 



^y (/(x, t)^) and the correlation length ^(t). Under the transformation L£, the correlation 
length scales like x such that ^ bS,. The surface width scales like / such that W — > b^W. 
These conditions are fulfilled if W and have a powerlaw dependence on time: 

m = Cot^ (1.10) 

W{t) = Wo t^ where /3 = ^. (l-U) 

Similarly, on mounded surfaces the average perimeter of the mounds increases oc t^^^ , while 
their height is oc t^. If the growth process selects one particular slope mo, the ratio of height 
and perimeter of the mounds is constant, which implies 1/z = /3 = a/z^a = l. 

Additionally, there are conditions for the functional form of the height-height correlation. 
Since dynamic scaling links properties of the surface at different times, we will use the notation 
G{l,t) for the height-height correlation at time t in the following. Since G{l,t) scales like the 
square of /, we have G{bl, bH) = b'^°'G{l, t). This is fulfilled if 



Gil,t) = f'^u -^\. (1.12) 




Since surface heights at distances greater than ^(t) are uncorrelated, G{l,t) saturates at large 
distance I. Consequently, u{x) ~ for large 

Similarly, the standard deviation w{L,t) of f{x,t) in a square of size L x L must fulfil 

w{L,t) = L-v (^-^y (1.13) 

For small t and large L where ^(4) <C L, the finiteness of L is irrelevant such that w{L,t) = 
W{t). This implies that v{x) ~ x^ for small x. On the other hand, for ^{t) ^ L, we expect 
w{L,t) to saturate such that v(x) ~ const, for large x. Then, w{L,t) assumes an asymptotic 
value 

Wl^te{L) OC L". (1.14) 

In Monte Carlo simulations, only systems of finite size can be considered. In a square 
system of size N x N, the deviations from the physically relevant case of an infinite system 
are small as long as ^(t) <^ N. As soon as the correlation length is comparable to N, finite 
size effects become relevant. In general, the properties of a finite system are similar to those 
of a section of an infinite system. Thus, the surface width of a finite system resembles the 
properties of w{N,t). 

So far, no assumptions on the morphology of the surface have been made. If we assume that 
the surface has self-affine properties on lengthscales smaller than ^(t), additional conclusions 



on the scaling behaviour can be drawn. Comparing equations |l.7| and |1.14, we find that the 
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exponent a and the Hurst exponent H must be identical. Additionally, from equation |1.5| we 
expect the height-height correlation G{l,t) to increase ~ P'^ for small /. Consequently, we 
have u{x) ~ const, for small \x\. In the literature, the dynamic scaling properties of self-affine 



surfaces are kown as Family-Vicsek scaling |34|. For historical reasons, deviations from this 



behaviour are denoted as anomalous scaling |6^, |9l[| . 

1.6 Atomic layer epitaxy 

Kinetic roughening and unstable growth make it difficult to grow perfectly flat layers by means 
of MBE. There is an alternative called atomic layer epitaxy (ALE) which, however, requires 
specific material properties |44]. ALE works for a binary compound AB which consists of 



alternating layers of A and B atoms which are parallel to the surface. The chemical bond 
between A and B atoms must be considerably stronger than A-A and B-B bonds. These 
conditions are fulfilled e.g. for (001) surfaces of II- VI semiconductors. 

In contrast to MBE of a compound material where the constituents are deposited simulta- 
neously, in ALE the surface is exposed to alternating pulses of pure A and B. Between the A 
and B phases, there is a dead time without any particle fiux. We start with the deposition 
of A on a flat 5-terminated substrate. Particle flux F and pulse time tp are chosen such 
that the amount of deposited material is sufficient for more than one layer of A. ^-atoms in 
the first layer are bound via the strong A-B bond. Atoms in the second or a higher layer 
have only weak A-A bonds. Temperature and dead time are chosen such that these surplus 
particles desorb, while atoms bound to the substrate remain on the surface. Since the rate at 
which atoms with A-B bonds desorb is much smaller than the desorption rate of atoms with 
A-A bonds only, no fine-tuning of the parameters F, tp and t^i is required for that purpose. 
Instead, we obtain the deposition of exactly one layer of ^ in a self-regulated manner. In the 
following 5-phase and the subsequent dead time, the same process repeats with the roles of 
A and B interchanged. At the end of the cycle, one monolayer of AB has been deposited. 
The whole procedure may be repeated an arbitrary number of times to obtain a film of the 
desired thickness. 

This behaviour might be complicated by the presence of surface reconstructions. In many 
compound semiconductors a surface terminated with a complete layer of one particle species 
is unstable. In particular, cation terminated surfaces of II- VI semiconductors like CdTe are 
characterized by vacancy structures with submonolayer coverage. Such effects may lead to 
ALE growth at rates 7^ IML/cycle. This feature has been used to investigate the stoichiom- 



etry of CdTe(OOl) and ZnTe(OOl) surfaces in experiments 35, 123]. In chapter 0, we will 



simulate ALE of reconstructed surfaces within the framework of our model of CdTe. 



1.7 Sublimation 

Occasionally, atoms at the surface of a crystal desorb. At a solid in thermal equilibrium 
with its vapour, this loss of material is compensated by the adsorption of an equal number 
of particles. In ultra high vacuum, this is not the case such that the crystal shrinks in the 
course of time. This non-equilibrium process is denoted as sublimation. 

The sublimation rate rgub is defined as the number of monolayers which desorb per unit 
time. If Tsub is measured as a function of T, one finds that in a wide temperature range 
rsub(^) is described well by an Arrhenius law (equation |1.2| ) with attempt frequency z/gub and 
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activation energy -Esub- Monte Carlo simulations |10C, 102| have shown, that the macroscopic 
parameters fgub and -Egub are not equal to the corresponding parameters of some rate limiting 
microscopic process like e.g. the desorption of a single adatom. Instead, they are functions 
of all model parameters. Consequently, to understand sublimation we have to consider the 
interplay of various physical processes on the surface similar to the investigation of crystal 
growth. 

At moderate temperature, the desorption rate of atoms bound to steps and islands can be 
neglected compared to that of single adatoms. Adatoms diffuse over a typical distance ^jes; 
the desorption length before they are desorbed. An adatom is created when an atom detaches 
from a step edge or when an atom jumps out of a closed layer. In general, the activation energy 
of the first process is smaller. At the latter process, a vacancy is created simultaneously. An 
atom at a neighbouring site in the same layer which diffuses into the vacancy leaves a new 
vacancy behind at the site where it came from. Conversely, one might say that the original 
vacancy has moved. Therefore, this process is called vacancy diffusion. 

Vacancies diffuse around until they are either filled by an adatom, are captured by a step 
edge or meet other vacancies and nucleate a vacancy island. This is a hole in the terminating 
layer of the crystal which exceeds some critical size such that it is unlikely to be closed again. 
Similar to the typical island distance li in submonolayer growth, in the sublimation of a flat 
surface there is a typical distance ly between neighbouring vacancy islands. Vacancy islands 
grow as they capture vacancies diffusing around and adatoms detach from their edges and 
desorb. Finally, they coalesce with their neighbours and one layer of the crystal is gone. This 
process is denoted as layer-by-layer sublimation. 

Similar to crystal growth, sublimation is influenced strongly by the presence of steps on 
a vicinal surface. If the terrace width Is is smaller than ly, the creation of vacancy islands 
is unlikely since vacancies are captured by the steps. Instead, sublimation proceeds in step- 
flow mode where the steps move back at a speed fgtep- At the desorption of one monolayer, 
the steps move over a distance Is such that Tgub = ^stcp/'s- Pimpinelli and Villain |^3| have 
investigated step-flow sublimation by means of BCF theory. 

If the spacing between steps is large compared to the desorption length, it is unlikely 
that an adatom which detaches from a step reaches a neighbouring step. Consequently, the 
interaction between steps is weak and Vstcp is independent of Ig. This case is denoted as free 
step- flow. In hindered step-flow, we have Is < Idcs- Then, there is a high probability that 
adatoms which detach from a step are captured by another step instead of being desorbed. 
Consequently, at small interstep spacing, Ugtep decreases if Is is reduced. 

Vicinal surfaces where the width of all terraces is equal sublimate either in layer-by-layer 
mode or in step-flow mode. However, techniques of surface preparation frequently yield a wide 
distribution of terrace widths. In this case, there might be a temperature range where ly is 
greater than the width of small terraces, but there are also large terraces which are wider than 
ly. Then, small terraces sublimate in step-flow mode. Simultaneously, there is layer-by-layer 
sublimation on large terraces. This effect has been observed experimentally on (001) surfaces 
of CdTe p|. 

In general, -Bgub depends on the mode of sublimation. Additionally, different surface 
reconstructions yield different values of -Bgub- This problem will be addressed in chapters ^ 
and |. 
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Chapter 2 

Flat (001) surfaces of 1 1- VI 
semiconductors: A lattice gas model 



Within the last years, potential technological applications of electronic devices based on II- VI 
semiconductors |116[] have inspired basic research concerning surfaces of these materials. In 
this context, various studies have addressed the properties of surface reconstructions. Exper- 

and how the 



imental studies have investigated which reconstructions are present |74, |75|, |1K 
reconstruction of the surface is influenced by parameters like temperature and particle flux in 
an MBE environment U7\, [I 



129|. The majority of this work has been devoted to CdTe and 



ZnSe, where a fairly complete qualitative overview over the phase diagram has been gained. 
An overview over the properties of CdTe can be found in [^]. On the other hand, there have 
been theoretical investigations of the reconstructions of CdTe [4C] and ZnSe 39, 79] using 
density functional theory. In these studies, knowledge about the chemical bonding of surface 
atoms and ground state energies of various reconstructions has been gained. 

In this and the following two chapters, we present models which aim at a theoretical under- 
standing of the reconstructed (001) surface of CdTe and ZnSe under typical MBE conditions. 
As explained in section 1.1, to study thermodynamic effects like phase transitions we need 
to consider systems with a large number of atoms. On the other hand, the numerical effort 
should not exceed the capabilities of today's computers. Therefore, efficient representations 
of atomic interactions are required. In many cases, two-dimensional lattice gases have been 
used successfully to model atoms adsorbed on a singular crystal surface or the terminating 
layer of such a crystal [jl^, |5^, |9^ . In spite of the conceptual simplicity of such models the 
interplay of attractive and repulsive short range interactions can result in highly nontrivial 
critical behaviour and complex phase diagrams. In this chapter, we will follow this approach 
to model flat (001) surfaces of CdTe and ZnSe in thermal equilibrium, our main focus being on 
CdTe. In the following two chapters, we will extend this model to consider three-dimensional 
crystals in nonequilibrium situations like growth and sublimation. 

The outline of this chapter is as follows: In section 2.1, we review the known facts about 
(001) surfaces of CdTe and ZnSe. In section 2^, we introduce a lattice gas model which 
considers the occupation of Cd sites and the dimerization of Te atoms and discuss its phase 
diagram. We conclude with a comparison of the features of our model with experimental 



results in section 2.3 
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2.1 Surface reconstructions of CdTe and ZnSe 



Both CdTe and ZnSe crystallize in the zinc-blende lattice. A sketch of the unit cell is shown 



in figure 2.1. This lattice is equivalent to a diamond lattice apart from the fact that the 



lattice sites are occupied alternately with the two elements. Each atom is surrounded by 
four nearest neighbours of the other species the nuclei of which lie on a tetrahedron. The 
zinc-blende lattice can be decomposed into two intersecting fee lattices the sites of which 
are occupied with particles of one species. Thus, it is composed of alternating half-layers of 
metal atoms (cations) and nonmetal atoms (anions) which are parallel to the (001) surface. 
Consequently, an ideal (001) surface would be terminated by a complete layer of one particle 
species. The positions of the atoms in one layer lie on a regular square lattice the axes of 
which are oriented in the [110] direction and the [110] direction. The lattice constant of this 
square lattice is a/\/2, where a is the lattice constant of the zinc-blende lattice itself. 




Figure 2.1: Unit cell of the zinc-blende lattice. 
Anions are shown in light grey, cations in dark 
grey. 

The orientation of the bonds of the atoms in one layer breaks the four-fold symmetry of 
the surface. Each surface atom is bound to two atoms of the opposite species in the half-layer 
below. In cations, these bonds point in the [110] direction. Conversely, surface anions have 
bonds in the [110] direction. Therefore, the symmetry of the (001) surface of the zinc-blende 
lattice is two- fold. 

The stability of the surface reconstructions can be understood from the electron counting 
rule [ 



pO| which was introduced in section |l.l| . Each metal atom provides two electrons. In 
each bond to a bulk atom, 1/2 electron is engaged. A nonmetal atom brings six electrons, 3/2 
of which are engaged in each bond to an atom in the bulk. The electron counting rule states 
that the remaining electrons must be sufficient to fill each anion dangling bond and every 
bond of a surface dimer with two electrons, whereas the cation bonds must remain empty. In 
particular, a surface terminated with a complete layer of undimerized atoms is unstable. 

Under vacuum, the (001) surface of CdTe is Cd terminated. The surface is characterized 
by vacancy structures where less than one half of the potential Cd sites in the top layer 
are occupied [27, 123|. At low temperature, one finds a c(2 x 2)cd reconstruction [^, |119|] , 
where Cd atoms and vacancies arrange in a checkerboard pattern (figure |2.2| a). Frequently, a 
contribution of a (2 x l)cd arrangement can be found. In the (2 x l)cd structure, the Cd atoms 
arrange in rows along the [110] direction which alternate with rows of vacancies (figure |2.2[ b). 
Density functional calculations [EH] have shown that the surface energies of the two competing 
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(a) (b) (c) (d) 



Figure 2.2: Panels (a),(b),(c): Sketches of the reconstructions of CdTe(OOl) which are dis- 
cussed in this work. Larger and smaller filled circles denote Cd atoms in the first and third 
layer, open circles represent Te atoms in the second layer. The grey rectangles show the sur- 
face unit cells. The [110] axis is aligned horizontally, (a), (b): Cd terminated reconstructions. 
Panel (a) shows the c(2 x 2)cd reconstruction, panel (b) the (2 x l)cd reconstruction. Panel 
(c) shows the (2 x 1)tg reconstrution of surfaces terminated with a complete monolayer of Te. 
Panel (d) shows the attractive couplings in our model. 



reconstructions are nearly degenerate and differ only by a small amount /S.E ~ 0.008 eV per 
(1 X 1) surface unit cellQ. At a temperature T ss 570 K, there is a phase transition | 119 | above 
which the (2 x l)cd arrangement of the Cd atoms dominates the surface. An analysis of high 



resolution low energy electron diffraction (HRLEED) peaks 75] has shown, that there is 
a high degree of disorder in this high temperature phase. One finds elongated domains with 
a large correlation length of ~ 375 A in the [110] direction and a domain width as small as 
22 A in the [110] direction. 

If CdTe is exposed to an external Cd flux in an MBE chamber, the c(2 x 2)cd reconstruction 
is stabilized at temperatures above the transition in vacuum. Under a Te flux, the surface is 
Te terminated with a (2 x l)Te reconstruction. At small Te fluxes, the surface is terminated by 
a complete monolayer of Te. The Te atoms on the surface form dimers which arrange in rows 
(figure |2.2| c). At high Te flux and low temperature, one additional Te atom is incorporated 
into each dimer, such that the Te coverage of the surface is 1.5. The symmetry of this 
reconstruction is still (2x1), since the Te trimers tend to arrange in rows. A schematical 
phase diagram of the surface following [^l| is shown in figure p.3|. 

Qualitatively, the properties of the reconstructions of ZnSe |2^, |12g| ] are quite similar to 
those of CdTe, where the Zn atoms are the counterparts of the Cd atoms and the Se atoms 
those of Te, respectively. There is one important exception: to date, no (2 x l)zn reconstructed 
surface has been found at high temperature. Density functional calculations [^, |7^ yield 
a higher energy difference IS.E ~ 0.03 eV per (1 x 1) surface unit cell between ideal c(2 x 2) and 
(2 X 1) reconstructed surfaces, which is approximately twice the value calculated for CdTe. 
As we will show, this greater energy difference might explain the different behaviour of CdTe 
and ZnSe. 



^In ^] a greater value of AE = 0.016 eV is given. However, AE is rather sensitive to some approximations. 
We use the unpublished result of a more precise calculation by S. Gundel. 
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Figure 2.3: Schematical sur- 
face phase diagram of CdTe 
in an MBE environment ac- 
cording to |21|. Shown is the 



predominant reconstruction as 
a function of temperature and 
particle flux. The properties 
of the surface in vacuum are 
shown in the middle. The up- 
per part displays the recon- 
struction under a Cd flux, the 
lower part shows the features 
of the surface under a Te flux. 



2.2 The lattice gas model 

The basic structure of our model of (001) surfaces is the same for different II- VI semicon- 
ductors. The differences between the materials are represented by the numerical values of 
parameters. In the following, for simplicity we will loosely speak of Cd and Te atoms instead 
of "cations" and "anions" without restricting ourselves to a modelling of CdTe only. 

We present a lattice gas model with effective interactions which include effects of surface 
strain and, at elevated temperature, lattice vibrations. Therefore, there is no simple mapping 
between ground state energies of the lattice gas model and surface energies determined from 
density functional theory. 

We model a flat (001) surface of CdTe, i.e. we neglect the influence of steps on the recon- 
struction. This is a reasonable approximation if the typical distance between steps is much 
greater than the size of the unit cell of the reconstruction. In thermal equilibrium, this is 
fulfilled at temperatures well below the roughening transition. In this case, in the absence of 
bulk vacancies or other defects the crystal is uniquely described by the state of the topmost 
monolayer which consists of one Cd half-layer and a Te half-layer. In our model, we consider 
such a monolayer with the Cd atoms on top. For simplicity, we assume the Te layer to be 
fully occupied. This is not a severe restriction, since the removal of Te atoms will uncover 
Cd atoms in the layer below. Within the limit of a model of a flat surface, we simply do not 
distinguish whether a Cd terminated surface has been created by removing a Te layer or by 
depositing additional Cd atoms onto an intact Te layer. 

We use cartesian coordinates where the x-axis points in the [110] direction, the y-axis 
points in the [110] direction and the z-axis points in the [001] direction. The origin of the co- 
ordinate system is at the centre of a Cd atom. Since we consider a flat surface, the z-coordinate 
of all atoms of one species is identical and will be omitted in the following. Measuring the 
lattice constant in appropriate units, the Cd atoms are at positions {x,y) with integer x,y. 
The Te atoms are at positions {x,y + 1/2), where x,y £ 7j. Tellurium dimers are created by 
the formation of a chemical bond between neighbouring Te atoms in y-direction ([110]). Since 
this is also the direction of the bonds of a surface Cd atom, a Te dimer is formed by a pair of 
atoms which might also be the binding partners of a Cd atom above them. This suggests a 
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simple lattice gas representation both of the occupation of Cd sites and Te dimerization: We 
consider a rectangular array x of integers {a^ijlf^^i- Xi^j = 1 represents a Cd atom at site 
which is bound to the Te atoms at sites (i, j — 1/2) and (z, j + 1/2). Xjj = 2 corresponds 
to a dimerization of Te atoms (i, j — 1/2) and + 1/2). Otherwise, Xjj- = 0. 

In principle, one might also consider the formation of Te trimers by introducing a fourth 
state Xjj = 3. This corresponds to a trimer which consists of the Te atoms at — 1/2) 
and + 1/2) and an additional Te atom at However, the formation of Te trimers 

plays a role only if the surface is exposed to a strong flux of pure Te at comparatively low 
temperature. Therefore, we have neglected this effect to reduce the number of parameters of 
our model and the numerical effort of its investigation. 

2.2.1 Interactions of atoms and dimers 

The detailed representation of the surface energies of a II- VI compound certainly would re- 
quire long-range interactions and terms which depend on the simultaneous occupation of three 
or more sites. However, it is plausible to assume that the dominant contribution to the surface 
energy stems from pairwise interactions of atoms at short distances. In this section, we intro- 
duce a Hamiltonian which considers pairwise interactions between Xij on nearest neighbour 
sites and diagonal neighbour sites. These reflect the known properties of the reconstructions 
of CdTe and ZnSe. 

Due to the electron counting rule, surfaces with Cd coverages greater than 1/2 are unstable 
while both a c(2 x 2)cd reconstruction and a (2 x l)cd reconstruction are permitted. This fea- 
ture can be captured by introducing a hardcore repulsion between Cd atoms on neighbouring 
sites in the y-direction. In the x-direction, an attractive interaction favours the occupation of 
nearest neighbour pairs the strength of which is denoted by Ex < 0. An attractive interaction 
Ed < between diagonal neighbours stabilizes the c(2 x 2)cd reconstruction. The parameters 
e-x and Ed are chosen such that the energy difference AE = \2eci — ex\/2 between these two 
reconstructions is small compared to the total surface energy per site. 

The electron counting rule favours Te dimerization, but forbids the formation of additional 
bonds of dimerized Te atoms. Therefore, we forbid the simultaneous occupation of neighbour 
sites in the y-direction with dimers (formation of chains of Te-Te-bonds) and Cd atoms next 
to a dimer. These rules permit both a (2 x l)xc reconstruction, where the dimers arrange in 
rows and a c(2 x 2)Te reconstruction, where they arrange in a checkerboard pattern. Density 
functional calculations [^] have shown that the surface energy of (2 x l)Te is significantly 
lower than that of a c(2 x 2)tc reconstruction. This may be understood from the more efficient 
relaxation of surface strain in (2 x 1)tc: Since the lattice is deformed in the same direction 
by dimers on neighbouring sites in x-direction, such an arrangement will be energetically 
favourable. We consider this fact by an attractive nearest neigbour interaction et between 
dimers in x-direction. Additionally, a dimer contributes a binding energy Sf, to the surface 
energy. Apart from the hardcore repulsion, we consider no interaction between Cd atoms and 
dimers. An overview over the couplings in our model can be found in figure [2.2|d. 

This yields the Hamiltonian 

L,N 

H = ^ ExCijCi^ij + EdCij (Cj+ij+i + Cj+ij-i) + Etdi^jdi^ij + Ehdij — HCij, (2.1) 

where we have introduced Cij = ^,1 and dij = ^,2- The boundary conditions are assumed 
to be periodic. Cij represents the occupation of lattice sites with Cd atoms: Cij = 1 at the 
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positions of Cd atoms and zero otherwise. Similarly, the positions of dimers are given by 
dij = 1. /Li is the effective Cadmium chemical potential, which includes the binding energy of 
surface Cd atoms to the substrate. The groundstate of the system at T = is determined by 
the chemical potential fj,. For fi > fiQ = 2ed — £t — ^b, the surface configuration with minimal 
energy is c(2 x 2)cd- For ji < hq the groundstate is (2 x l)Te- 

In the following, we measure energy in dimensionless units which have been adjusted such 
that £d = —1. Additionally, we set Boltzmann's constant k = 1 such that temperature is 
measured in units of \ed\- 

2.2.2 Characterization of the surface configuration 

To characterize the surface reconstruction quantitatively, we evaluate the mean Cd coverage 
PCd = and the correlations 

Ccd = ^ (ci,i (ci+i,i+i + cmj-i))y (2.2) 
C^d = (2.3) 

Cqj^ measures the probability to find two diagonal neighbour sites which are simultaneously 
occupied by Cd atoms. This is a measure for the fraction of the surface which is covered 
by regions with a local c(2 x 2)cd order. Similarly, C^^ measures the contribution of locally 
(2 X l)cd reconstructed regions in the system. 

The long range order of the Cd atoms is measured by the order parameters 

L,N 
^ L,N 

Mga'^ = (2-5) 

^Cd"^^ is the staggered magnetization of a system of Ising variables {sijjf where Sjj = 1 
if Cjj = 1 and —1 otherwise. Large absolute values indicate a long range order of the c(2 x 2)cd 

reconstruction. Its counterpart M^^^^ measures the long range order of (2 x l)cd- 

Further, we introduce the mean Tellurium dimerization p£) := {dij)- ■ and the dimer corre- 
lation Cf) := (dijdi+ij)- J, which characterizes the number of dimers which are incorporated 
in locally (2 x 1) ordered areas. Their long range order is characterized by 



M 



L,N 

(2x1) 

^ LN ^ 



TT? yZ ^id COS (vrj) . (2.6) 



2.2.3 Methods of investigation 

We have investigated this model by means of the transfer matrix method and Monte Carlo 
simulations. An introduction to the transfer matrix technique is presented in appendix 
This method allows for a numerical calculation of the free energy of lattice systems with short- 
range interactions. In general, the system size must be finite in all directions but one. In our 
investigations, we choose the y-direction as the infinite direction. Since the computational 
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effort increases exponentially fast with the system size L in the x-direction, only comparatively 
small L are feasible. However, this is sufficient to investigate the phase transitions of the 
system in the limit L ^ oo via finite size scaling techniques. To this end, we have applied the 
method discussed in section A.2.2 to determine the loci of first order phase transitions and 
to compute coverage discontinuities. Continuous transitions have been investigated by means 
of phenomenological renormalization group theory, which will be introduced in section |A.2.1 



Coverage and correlations can be obtained from proper derivatives of the free energy or directly 
from the relevant eigenvector of the transfer matrix. We have refrained from calculating order 
parameters by means of the transfer matrix. This would require the introduction of staggered 
fields in the Hamiltonian 2.1. A system with this modified energy function is numerically 
intractable for reasonable values of L. 

Additionally, we have performed Monte Carlo simulations of our model. To obtain a rea- 
sonably fast equilibration, we have applied continuous time algorithms which will be discussed 
in more detail in appendix We have simulated both the grand-canonical ensemble where 
Pcd is controlled by fi and the canonical ensemble where the number of Cd atoms in the sys- 
tem is fixed. In general, the results of the transfer matrix calculations and the Monte Carlo 
simulations are in good agreement (figures 2.4 and |2.6|). 



2.2.4 Simplified model without Te dimerization 

In this section we present a simplified version of our model where we do not consider the 
dimerization of Te atoms, i.e. the possible values of the OC 2 J 3)1* G constrained to the values 
and 1. Then, the state of the system is uniquely described by {cij}^^^^. The Hamiltonian of 
the simplified model is 

L,N 

Hs = ^ ExCijCi+i^j + EdCij (Cj+ij+i + Ci+ij_i) — HCij. (2.7) 

For fi > —2, the groundstate is a c(2 x 2)cd reconstruction, whereas for fi < —2, the 
system is empty at zero temperature. Figure |2.4| shows the temperature dependent behaviour 
of the correlations and order parameters of this model with the parameter set = — 1, 
Ex = —1.96 in the grand canonical ensemble at a chemical potential fj, = —1.96. In the Monte 
Carlo simulations, the system was initialized with a perfect c(2 x 2)c(j reconstruction at all 

investigated temperatures. At low temperature, we obtain values of M^^ and C^^ close to 

f 2 X 1) 

0.5. is small and is zero apart from small fluctuations which are finite size effects. 

This is the signature of a long-range ordered c(2 x 2)cd reconstructed phase. At T ~ 0.3, the 
system loses its long-range order such that at high temperature we have M^^j^^^ ~ M^j^^^ ~ 0. 
Simultaneously, the Cd coverage jumps to a smaller value. These features indicate a first order 
transition from the c(2 x 2)cd phase at low T to a disordered regime at high temperature. 

At the phase transition, Cq^ decreases rapidly. Conversely, displays a sudden increase 
such that C^^ > in the high temperature regime. Consequently, in the vicinity of the 
phase transition not only the long range order, but also the local arrangement of the Cd atoms 
changes. In the disordered phase they prefer to arrange in rows which are characteristic of 
the (2 X l)cd reconstruction. 

In figure \2,.5[ we show phase diagrams of the simplified model at Ex = —1.6 (figure ^^a) and 
Ex = —1.9 (figure ^b) in the pcd-T plane. In addition to the transfer matrix extrapolation, 
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Figure 2.4: Temperature dependent behaviour of the simplified model at = —1, Sx = fJ- = 
— 1.96. Panel (a): Cd coverage pcd {+), correlations C^^ (x), C^^ (^) and mean absolute 
of the order parameter M^^'^^ (©). The symbols show results of Monte Carlo simulations 
of a 128 X 128 system. lO'^ • LN events have been performed both for equilibration and for 
measurement. The lines show results of a transfer matrix calculation at a strip width L = 18. 
Panels (b), (c) show sections of 16\/2 x 16\/2 zinc-blende lattice constants of the surface at 
the end of the simulation run at T = 0.29 (b) and T = 0.33 (c). 




Figure 2.5: Phase diagrams of the simplified model without Te dimerization in the pcd-T 
plane at Sx = —1.6 (panel (a)) and £x = —1.9 (panel (b)). In both figures we have = — 1. 
The solid lines show the phase boundaries which have been extrapolated from transfer matrix 
calculations with strip widths L = 14, 16 and 18. The crosses show results of Monte Carlo 
simulations of a 128 x 128 system at fixed pcd- Phase (I) is homogeneously ordered with a 
c(2 X 2)cd reconstruction whereas phase (II) is disorderd. In region (III), phases (I) and (II) 
coexist. At the dashed line, we have C^^ = C^^. The symbols (x) show the respective results 
of the Monte Carlo simulations. Statistical errors of the simulation data are on the order of 
0.05. 
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we obtain a rough estimate of the phase diagram from canonical Monte Carlo simulations at 
constant Cd coveragen. For this purpose, a non-local algorithm which exchanges empty and 



occupied sites according to Kawasaki rates 126 1 is applied. A rapid decrease of M^^""^^ 



with T indicates the transition into the disordered phase. A search for a pronounced maximum 
in the fluctuation of the order parameters yields results which are identical within errorbars. 



The symbols in figure 2.5 show the results for N = L = 128 which are in reasonable agreement 
with the transfer matrix prediction. 

At low temperature and high pcd, the system is in the homogeneously ordered phase (I). 
The disordered phase (II) is found at extremely low Cd coverage. In the regime (III) at inter- 
mediate Cd coverages, the phases (I) and (II) coexist, pcd in both phases is given by the right 
and the left boundary of region (III), respectively. The coexistence regime (III) disappears 
at a tricritical point. At high temperature the system is either homogeneously ordered or 
disordered, depending on coverage. The transition between these regimes is continuous in 

(2x2) 

terms of and should belong to the Ising universality class. Our transfer matrix calcu- 

lations yield a temperature exponent^ = 1 which is consistent with this picture. However, 
a calculation of the exponent yh is not possible without the introduction of staggered fields. 
In the limit T ^ oo, the line of the transition between the phases (I) and (II) approaches 
the line pcd = 1/2- This can be understood from the fact that at very high temperature the 
infinite repulsion between nearest neighbours in y-direction is the only relevant interaction. 
Then, colums of lattice sites aligned in y-direction decouple such that the system is disordered 



at arbitrary pcd- This is in contrast to isotropic hard square lattice gases ||lj, |5l|] where an 
extended regime (I) persists for arbitrary temperature. 

To characterize the short range order of the Cd atoms, we have determined the line where 
^Cd ~ ^Cd- Since this is not a line of phase transition, it is not possible to apply finite 
size scaling theory to extrapolate its position to the limit of an infinite system. However, the 
results of transfer matrix calculations for various L and the canonical Monte Carlo simulations 
at = L = 128 agree well such that finite size effects are probably weak. At the right of the 



dashed line in figure 2^, the c(2 x 2)cd ordering is prevalent. At smaller pcd^, the Cd atoms 



arrange preferentially in a (2 x l)cd order. 



Figure |2.4| demonstrates the crucial role that the difference between the surface energies of 
perfect c(2 x 2)cd and (2 x 1)^^ reconstructions plays for the phase diagram. With increasing 
/\E, the tricritical point shifts to smaller coverage and higher temperature. Even more so 
does the line which separates c(2 x 2)cd from (2 x l)cd prevalence. 

2.2.5 Results of the model with Te dimerization 

To get insight into the typical behaviour of the model, we will first present a detailed investi- 
gation using the parameter set = = — 1, e^; = ej = — 1.9. The choice Ex = —1.9 reflects 
the fact that the energy difference between c(2 x 2)cd and (2 x l)cd is small. The choice 
Efe = — 1 yields low concentrations of Te atoms which are neither dimerized nor bound to Cd 
atoms. This is consistent with experimental results. For simplicity, we have chosen £t = £x 
as a typical example of the case where both couplings are of the same order of magnitude. 
However, since there is no next-nearest neighbour interaction between dimers, there is a sig- 
nificant difference in the surface energies of (2 x l)Te and c(2 x 2)tc, which is consistent with 

^These simulations have been performed by T. Volkmann. See ]l26{ for additional information. 
^The exponents yT and yn are defined in appendix A. 2.1 
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the results of |0|]. 

Below, we discuss the influence of the numerical values of the parameters on the phase 
diagram. We will show, that a smaller value of ej affects the phase diagram qualitatively. In 
particular, this concerns the high temperature phases. 

Grand-canonical ensemble 

With the above parameter set, the Cd and the Te terminated groundstates are separated 



by Mo = 0.9. In figure 2^ the temperature dependent behaviour of the system at constant 
chemical potential is shown for fJ, = 1 (figure p.6| a,b) and /i = 0.8 (figure ^I^c,d), which are 
examples for both cases. At ^ = 1, the ordered c(2 x 2)cd phase at low T manifests itself in a 

, (2x2) 

high Cd coverage pcd find values of the correlation C^^ and the order parameter ' close 



to 0.5. In figure |2.7| d a surface snapshot in this phase at T = 0.71 is shown. As expected, 
the Cd atoms arrange preferentially in a checkerboard configuration. At a temperature T = 
0.83, there is a first order phase transition to a disordered phase. This is indicated by a 

(2x2) r1 

decrease of the order parameter and the correlation C^^. Simultaneously, C^^ starts 

to increase. At T = 1, both lines cross such that the high temperature behaviour of the system 
is dominated by a local (2 x 1) ordering of the Cd atoms. The simulation data shown in figure 



2.6|a,b have been taken from two independent simulation runs. The system was initialized 
with a perfect c(2 x 2)cd configuration at the lowest temperature investigated. Then, as 
successively higher temperatures were imposed, the surface configuration was kept as initial 
state of the next simulation. This reflects in a small hysteresis effect in due to the first 

order nature of the phase transition. Qualitatively, the behaviour of the Cd atoms is quite 



similar to that found in the simplified model (figure 2.4). 



At /U = 0.8, we obtain a completely different behaviour. Since the groundstate is a (2 x l)Te 
reconstruction, at small T we measure low Cd coverages and values of C|j and close to 

0.5. The most frequent thermal excitations are Cd adatoms, the density pcd of which increases 
with T. Figure |2.7| c shows a surface snapshot at T = 0.93. The Cd atoms preferentially 
arrange in rows such that Cq^ ^ ^cd- These rows adapt to the structure of the (2 x l)Te 
reconstruction. This yields nonzero values of the order parameter M^^^\ which indicate a 
global ordering of Cd atoms in a (2 x 1) arrangement. However, the interactions between 
the Cd atoms themselves are insufficient to stabilize this global order. Instead, it is purely 
induced by the environment of the Te dimers. 

At T = 1.23, there is a first order phase transition above which the system is in a disordered 
phase similar to that found at // = 1. Remarkably, we observe high pcd ~ 0.2 at temperatures 
slightly below the phase transition. At the phase transition, pcd jumps to an even higher 
value. Above the transition it decreases slightly with T. 

Phase diagram 



Figure 2.7a,b shows the phase diagram of the model, which has been extrapolated from 



transfer matrix calculations with strip widths L of 6, 8 and 10. In figure 2.7b, the lines of 
phase transitions in the p-T plane have been plotted. At low temperature, the system is 
either in an ordered (2 x l)Te (1) or an ordered c(2 x 2)cd phase (2). The line of phase 
transition between these phases starts at zero temperature and p = pQ where the energies of 
both reconstructions are degenerate. For < T < Tt = 0.84, it remains at the same chemical 
potential pQ, apart from small numerical uncertainties of extrapolation. In consequence, a 
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Figure 2.6: Simulations at constant chemical potential. Panels (a), (b) show the behaviour 
of the system with the couplings £d = £b = — Ij = £t = —1-9 and a cadmium chemical 
potential ^ = 1. With these parameters, the ground state of the system is c(2 x 2)cd- Panel 

(a) shows the Cd coverage pcd (+) and the correlations C^^ (x), (^) and (□). Panel 

(b) shows the mean absolute of the order parameters M^^"^^ (0), M^^^^ (■) and M^^^^ (•). 
The symbols show data from a simulation run at a system size L = N = 64. 4 • 10*^ • LN 
events have been performed both for equilibration and for measurement. The lines have been 
obtained by means of a transfer matrix calculation at a strip width L = 10. The data shown 
in panels (c), (d) have been obtained at /i = 0.8, where the ground state is (2 x l)Te- All 
other parameters are identical to those used in (a), (b). The meaning of the symbols is the 
same as in panels (a) and (b). 
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Figure 2.7: Phase diagram witii parameters ea = £b = — !> £x = = —1-9. Tlie solid 
lines in panel (a) show the lines of phase transitions as functions of pcd and T, in panel (b) 
they are shown in the fi-T plane. Note the offset in the temperature axes. In region (1), the 
system is in a homogeneously ordered phase with a (2 x 1)tc reconstruction, in region (2) it 
is homogeneously ordered and c(2 x 2)cd reconstructed. Region (3) is the disordered phase. 
On the left side of the dashed line, the Cd atoms show preferentially a local (2 x 1) ordering 
(^Cd ^ ^Cd)^ while on its right side a c(2 x 2) ordering dominates. Due to the coverage 
discontinuity between (1) and (2), (3), there is a coexistence regime where regions with high 
and low pcd coexist (4). The symbols show lines of constant /i = 0.8 (^) and /i = 1 (□) (same 
data as in figure |2.6|a,c). Panels (c)-(f) show typical surface snapshots, (c)-(e) correspond 
to grand-canonical simulations. The snapshots show surface configurations after 8 • 10^ • LN 
events, (c): T = 0.93, /x = 0.8 (phase 1), (d): T = 0.71, p = 1 (phase 2), (e): T = 1.14, 
/i = 1 (phase 3). (f) displays a surface configuration after 2 • 10"^ • LN events in a canonical 
simulation at T = 0.73, pcd = 0.35 (Region 4, coexistence of phases 1 and 2). All snapshots 
show sections of x 16\/2 zinc-blende lattice constants of systems of size L = N = 64. 
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phase transition between a c(2 x 2)cd and a (2 x l)Te reconstruction cannot be observed 
at constant chemical potential. The disordered phase (3) exists for T > Tt- At the point 
(T = Tt,ii = /ig) five phases coexist: the disordered phase, and c(2 x 2)cd and (2 x l)Te 
reconstructed phases in two sublattices, corresponding to positive and negative values of the 
order parameters M^^"^^ and M^^^\ respectively. The (2 x 1)to reconstructed phase (1) 
exists only at temperatures below a critical temperature ~ \et\. At T}, the line of the 
phase transition to the disordered phase (3) diverges to /i = — oo. On the contrary, the 
c(2 X 2)cd reconstructed phase (2) may exist at arbitrary temperature, if the Cd chemical 
potential is large enough. The dashed line in figure |2.7| b shows the chemical potential at 
which the correlations Cq^ and Cq^ in the disordered phase are equal. For smaller /x, the 
local ordering of Cd atoms is dominated by a (2 x 1) arrangement, while for larger /i they 
prefer a local c(2 x 2) ordering. 

The phase transition from phase (2) to phase (1) occurs at a temperature independent 
chemical potential /iq- At the transition between the phases (2) and (3), pcd varies contin- 
uously. Therefore, there is no coverage discontinuity if temperature is varied at a constant 
chemical potential /j. > /j-q where the groundstate is a c(2 x 2)cd reconstruction. 

Figure |2.7| a shows the phase diagram in the pcd-T plane. There is a coverage discontinuity 
at the phase transition between the ordered phases (1) and (2) and at the transition between 
phases (1) and (3). These discontinuities yield a coexistence regime (4). Here, for T < 
Tt, c(2 X 2)cd and (2 x 1)tc reconstructed phases coexist, while for T > Tt the (2 x l)Te 
reconstructed phase coexists with the disordered phase. Figure 2^f shows a surface snapshot 
from a simulation which was performed at a constant Cd density pcd = 0.35 and a temperature 
T = 0.73. This is a typical surface configuration in the regime where the ordered phases (1) 
and (2) coexist. The system is separated in two phases, one with a high pcd and a c(2 x 2)cd 
reconstruction, and another one with a (2 x 1)tc reconstruction and a low concentration of Cd 
adatoms. The local values of pcd in both regions are given by the left and the right boundary 
of the coexistence regime. At low temperature one obtains pcd ~ 0.5 in the c(2 x 2)cd phase 
and pcd ~ in the (2 x l)xe phase. At temperatures T ^ 0.6, pcd in the (2 x l)Te phase 
increases strongly with T and obtains its maximal value 0.23 at T = 0.98. At even higher 
temperature it decreases with T and becomes zero at T^ . On the contrary, pcd in the c(2 x 2)cd 
phase remains high for T < Tt- At T > Tt, the Cd-rich phase is disordered. Then, pc^ at 
the right boundary of the coexistence regime decreases with T. At T^, it becomes zero and 
the coexistence regime disappears. The dashed line in figure 2.7a marks the values of pcd at 
which C^^ = C'q^. For smaller coverages, > C'^^ such that the local ordering of the Cd 
atoms is dominated by a (2 x 1) arrangement, while for greater coverages C^^ < C^^. The 
values of pcd which have been measured in the simulations shown in figure 2^ are plotted in 
the phase diagram as examples of lines of constant chemical potential. The loci of the phase 
transitions and the point where = C'^^ (at p = 1) are in good agreement with the results 
of the transfer matrix extrapolation. 



Canonical ensemble 

The temperature dependent behaviour of the model at a constant Cd adatom density pcd = 
0.35 is shown in figure |2.8| . This is an example which shows the typical behaviour of the 
model under the conditions of phase separation. At low temperature, most of the Cd atoms 
are concentrated in a Cd rich phase with a c(2 x 2)cd reconstruction. The remaining area of 
the system is covered with a (2 x l)Te reconstructed phase. Since both phases are long-range 
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Figure 2.8: Results of a simulation of a 64 x 64 system at conserved pcd = 0.35 and couplings 
£d = Efe = — 1, = = —1-9. 10^ • LN events have been performed both for equilibration 
and for data sampling. Panel (a) shows the correlations C^^ (x), C^^ (^) and (□). Panel 

(b) shows the mean absolute of the order parameters M^,^ (0), M^,^ (■) and (•). 
The inset shows standard deviations of order parameters. The meaning of the symbols is the 
same as in the large picture. 



f2x2) f2xl) 

ordered, we obtain large values of the order parameters , Mj^ and the corresponding 

correlations. The fraction of Cd atoms which occupy sites in the Te rich phase yields nonzero 
values of the order parameter M^^^\ As temperature increases, more and more Cd atoms 

(2x1) 

pass into the Te rich phase. This yields an increase in C^^ and and a decrease of 

C(^^ and due to the (2 x 1) ordering of the Cd atoms in the Te rich phase. At 

(2x2) 

the temperature Tt, the Cd rich phase undergoes the order-disorder transition where 
drops to zero. At this phase transition, the fraction of Cd atoms which are incorporated in 
locally (2 x 1) ordered configurations in the Cd rich phase increases. Thus, there are two 
independent effects which lead to a dominance of C^^ over C^^ at higher temperature: The 
specific shape of the left boundary of the coexistence regime and the order-disorder transition 
of the Cd rich phase. At a temperature T^, the system is leaving the coexistence region into 
phase (3). At this phase transition, the separation between a Te rich and a Cd rich phase 
disappears, such that at high temperature the system is in a homogeneous, disordered state. 
The order parameters M^^^^ and M^^^^ become zero and the local correlation between 

(2 X 1) 

dimers, Cj-, , decreases. In the simulations we have determined the critical temperatures of 
the phase transitions from the standard deviations of the order parameters M^^'^^ and M^^^^ . 

(2x2) (2x1) 

M^^ is peaked at Tt, where the long range order of the Cd atoms is lost, while is 
peaked at T^, where the system leaves the coexistence regime and the dimers lose their long 
range order. We obtain Tt = 0.79 it 0.2 and T^ = 1.0 it 0.2, where the uncertainty is due 
to the temperature spacing between the single simulations. These values are systematically 
lower than the theoretical results of the transfer matrix extrapolation (Tt = 0.84, T^ = 1.11). 
However, the theoretical results are valid in the limit of an infinite system size, where the free 
energy of the phase boundary can be neglected compared to that of the bulk of the phases. 
We have verified that the systematic deviation between theory and simulations decreases with 
the system size. 
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Figure 2.9: Phase diagram with the parameter set Ed = £b = — li = —1-6, et = —1.2. 
Compared to the phase diagram shown in figure the energy difference between the c(2 x 
2)cd and the (2 x l)cd reconstruction is greater and the interaction between dimers is weaker. 
With these parameters, the transition between phases (1) and (3) is continuous. Panel (a) 
shows the phase diagram in the pcd-T plane, in panel (b) it is shown in the fi-T plane. 



The influence of the parameter set on the phase diagram 

Finally, we discuss the influence of the parameter set on the phase diagram. The main effect 
of a variation of the binding energy of the Te dimers Si, is to shift the chemical potential at 
which the phase transitions occur. The influence on the shape of the phase diagram is small 
as long as £b is sufficiently large such that the density of Te atoms which are neither dimerized 
nor bound to Cd atoms is small. Due to the electron counting rule, a state where a large 
fraction of the Te atoms remains unbound should be irrelevant for CdTe. 

A smaller energy difference AE between a perfect c(2 x 2)Qij^ and a perfect (2 x l)cd 
reconstruction increases the tendency of Cd atoms to arrange in a (2 x 1) order. If = —1.95 
and all other parameters are identical to our standard parameter set, the dashed line where 
^Cd ~ ^Cd shifted to a higher pcd ~ 0.45. Additionally, we obtain a higher value of 0.28 
for the maximal Cd coverage in the Te rich phase, which is achieved at a lower temperature 
T = 0.9. The temperature Tj = 0.7 above which the disordered phase is stable is also slightly 
lower. In general, a greater value of \ex\ increases both the tendency of the Cd atoms to 
arrange in a (2 x 1) order in the disordered phase and the concentration of Cd atoms in 
the (2 X l)Te reconstructed phase. The temperature of the order-disorder transition in the 
Cd rich phase is lowered. Conversely, a greater energy difference between perfect c(2 x 2)cd 
and (2 x l)cd reconstructed phases alters the temperature of this transition and leads to a 
preferential c(2 x 2) arrangement of the Cd atoms in the disordered phase. 

The main effect of a variation of the interaction £t between Te dimers is to change the 
properties of the transition between the phases (1) and (3). As long as \£t\ is sufficently 



high, the phase diagram is qualitatively similar to that shown in figure 2.7. In this case, 
the main effect of a variation of £t is to shift the temperature ~ |ej| where the (2 x l)Te 
phase vanishes. However, at low \£t\ this transition becomes a continuous phase transition 



without any coverage discontinuity. As an example, in figure |2.7| a phase diagram with et = 
— 1.2, £ii = £d = —1 and £x = —1.6 is shown. The coexistence regime (4) vanishes at 
Tt = 0.9 such that there is no phase separation between the long-range ordered (2 x l)Te phase 
and the disordered phase. This parameter set yields a comparatively high AE. Therefore, 
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the line where = Cq^ is at relatively low Cd coverages. In contrast to the situation 
at large \et\, the chemical potential at which the transition from phase (1) to phase (2) 
occurs is not independent of temperature. Instead, the transition is at slightly greater n at 
higher temperature. In consequence, there is a small range 0.2 < fi < 0.27 where there is a 
phase transition from phase (2) to phase (1) if T is increased at constant chemical potential. 
Qualitatively, the phase diagram at small \£t\ is reminiscent of that of the simplified model as 
far as the arrangement of Cd atoms is concerned. 



2.3 Comparison with experimental results 



The results presented in section 2.2.5| suggest an interpretation of the experimentally observed 



crossover from a c(2 x 2)cd reconstruction to a (2 x l)cd reconstruction as an accompanying 
effect of an order-disorder phase transition. At low temperature, there is a long-range ordered 
Cd-rich phase with a c(2 x 2)cd reconstruction. At a critical temperature, the Cd atoms lose 
their long-range order and arrange preferentially in a (2 x 1) pattern. This picture is consistent 
with the experimental observation of small domains in the (2 x l)cd reconstruction |74, 75] 
which indicate a high degree of disorder. 

Strictly speaking, a CdTe surface under vacuum is not in thermal equilibrium. At the 
temperature of the c(2 x 2)cci - (2 x l)cd transition, sublimation plays an important role. 
As we will show in chapter |3|, important features of the simplified model which neglects Te 
dimerization are preserved under the conditions of step-flow sublimation. This is the dominant 
sublimation mechanism for CdTe (001) ||7^. However, the Cd coverage pcd is determined by 
the sublimation process. Therefore, within the limit of an equilibrium model it is not possible 
to calculate the path in the phase diagram that CdTe follows. 

Auger measurements ]7^ , |75| ] yield pcd ~ 0.35 at the transition temperature in vacuum. 
This suggests that the system is in the coexistence regime for a wide range of temperatures. 
Then, the behaviour of CdTe in vacuum should be similar to our results at constant pcd = 0.35 
(figure |2.8| ) . Electron diffraction techniques investigate large regions on the surface and hence 
yield averages over all coexisting phases. Consequently, it should be reasonable to compare 
results of these experiments with quantities which are averaged over the whole system in 



our canonical simulations. As discussed in section 2.2.5, there are two independent effects 
which lead to a clear dominance of the arrangement of Cd atoms in rows: the order-disorder 
transition in the Cd-rich phase and the fact that an increasing fraction of the Cd atoms is 
dissolved in the Te-rich phase. Electrons are diffracted both from the Cd atoms and the 
Te dimers which have a (2 x 1) order. The superposition of both effects should yield the 
pronounced (2 x 1) diffraction peaks which have been observed in [74, ITT 



In the c(2 x 2)q^ reconstructed phase, at temperatures well below the phase transition one 
frequently finds collective thermal excitations of adatoms, where a row of several Cd atoms is 
shifted by one lattice constant in the y-direction. Due to the repulsion between Cd atoms in 
this direction, this is possible only if one Cd atom per excitation is missing. An example is 



shown in figure p.7|d. This effect has been observed experimentally by Seehover et. al. [108] 
at room temperature using STM microscopy. There is a striking similarity between figure 3 
in l|l^ and figure l^d. 



The phase diagram of our model explains the properties of the reconstructions of the CdTe 
(001) surface under an external particle flux. If the Cd coverage of the surface is increased 
by deposition of Cd, the state of the system moves into regions of the phase diagram where 
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the Cd atoms arrange preferentially in a c(2 x 2) pattern. Depending on temperature, this is 
either the ordered c(2 x 2)cd phase (2) or the region of the disorded phase (3) on the right 
side of the dashed line where C^^ > C^^. Indeed, experiments |119[ have shown that an 
external Cd flux restores the c(2 x 2)cd reconstruction at high temperatures where a (2 x 1) 
order is found under vacuum conditions. Clearly, a strong Te flux induces the formation of a 
long-range ordered (2 x l)Te phase at temperatures below T]:. 

On ZnSe, a Zn terminated (2 x 1) reconstruction has not been observed yet. Our model 
offers an explanation of this fact which is consistent with the results of density functional 



theory ||3^, 39, 40, 7S]. These calculations show, that for ZnSe the difference between the 



surface energies of perfect cation terminated c(2 x 2) and (2x1) reconstructions is significantly 
greater than for CdTe (0.03 eV per (1 x 1) surface unit cell versus 0.008 eV). In our model, 
this greater energy difference corresponds to a smaller value of \£x\ which shifts the line where 
^Cd ~ ^Cd smaller coverages. Consequently, a c(2 x 2)zn arrangement dominates in a 
much wider range of coverages. 

Wolfframm et. al. [|l29f| have measured the locus of the transition beween the c(2 x 2)zn 
reconstruction and the (2 x l)ge reconstruction as a function of temperature and the com- 
position of an external particle flux by means of reflection high energy electron diffraction 
(RHEED). They find that at higher temperature a greater Se flux is required to obtain a 
(2 X 1) diffraction pattern. At temperatures above 450°C no (2 x l)se reconstruction could 
be observed even under extremely Se-rich conditions. This is reminiscent of our observation 
that the anion-rich phase (1) vanishes at a temperature T^. 

Unfortunately, the available experimental data are insufficient for a systematical fit of the 
model parameters. However, some rough estimates show that at least the orders of magnitude 
are reasonable. As discussed above, the c(2 x 2)cd -(2 x l)cd transition of CdTe in vacuum 
should be at a temperature close to Tt. Identifying this with the experimental value of 570 K, 
we obtain the value of our energy unit \eii\ ~ 0.06 eV in physical units. This yields a value 
/S.E ~ 0.003 eV for the difference in the surface energies of c(2 x 2)cd and (2 x l)cdj which 
is a bit less than 1/2 of the value of 0.008 eV which has been obtained by means of DFT 
calculations. This shows at least that the qualitative agreement between experiments and our 
model has been obtained in a physically reasonable region of the parameter space. 

In our model, phase (1) vanishes at a temperature w \et\. Identifying this with the value 
of 450°C which has been measured in experiments on ZnSe, we obtain that |e(| ~ 0.06 eV. 
This is the same order of magnitude as our estimate of \£d\ in CdTe. 

These considerations suggest that it should be possible to obtain quantitative agreement 
between experiments and our model both for CdTe and ZnSe with values of the model pa- 
rameters on the order of magnitude of a few ten meV. 
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Chapter 3 

Reconstructions on non-equilibrium 
crystal surfaces 



In section ^, we have presented a lattice gas model which displays an order-disorder phase 
transition the properties of which are reminiscent of the c(2 x 2)q(^-{2 x l)cd reordering of the 
Cd atoms at the CdTe(OOl) surface. However, a crystal surface in an MBE chamber is not 
in thermal equilibrium. In the presence of a particle flux, the crystal grows. Under vacuum, 
there is sublimation (section |1.7D . 

In the literature, a few models of crystal growth have been presented which take into 
account the effects of surface reconstructions. On the one hand, there have been models which 
aim at a detailed understanding of the features of specific materials, in particular Si 94] 



and GaAs |45, 56]. On the other hand. Chin and den Nijs ]2C] have investigated a simple 
model to address the fundamental question whether there might be a reconstruction order on 
growing surfaces. They find, that a long range order in the strict sense which is present in 
thermal equilibrium ]|3^ is absent during growth. However, the order of the reconstruction is 
preserved on small and intermediate lengthscales. 

To the best of our knowledge, phase transitions between different reconstructions have not 
been considered in growth models, so far. Consequently, a theoretical understanding of the 
interplay of this phenomenon with the dynamics of growing or sublimating surfaces has not yet 
been achieved. We address this important problem in the context of a model of a compound 
material that displays a competition of different vacancy structures in the terminating layer 



similar to the lattice gas model which was introduced in section 2.2.4. In particular, we 



address the following questions: (1) It has been claimed ]74, 75] that the surface layer is in 
effective thermal equilibrium under the non-equilibrium conditions of sublimation. Is this 
justifed? If not, which features of an equilibrium phase transition are preserved? (2) What is 
the effect of an external flux of particles of one species on the surface reconstruction? These 
simulations model situations which are given in experiments where phase diagrams like that 
shown in figure |2.3| are determined. 



3.1 The model 



Our model is an extension of the simplified planar lattice gas model introduced in section 2.2.4 
to a model of a three-dimensional crystal. In this chapter, we aim at an understanding of 
fundamental properties of nonequilibrium crystal surfaces rather than a realistic modelling of 
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CdTe. To this end, we consider a lattice gas model of a compound "^5" with a cubic lattice 
and a comparatively simple potential energy surface. The A (B) atoms correspond to Cd (Te). 
We impose the constraint that there are no overhangs and dislocations. This solid-on-solid 
condition makes is possible to describe the configuration of the crystal uniquely by a two- 
dimensional array h of integers {/ijj jf^-^^ which denote the height of the column of atoms at 
x-coordinate i and y-coordinate j. Thus, we achieve an effectively two-dimensional description 
of a three-dimensional crystal which reduces the complexity of our model considerably. 




Figure 3.1: Sketch of a surface showing the inter- 
actions between the particles. A particles are shown 
as dark cubes, B particles in light grey. The x-axis 
points from left to right, the y-axis from front to back. 



To model the layered structure of the CdTe crystal, we assign the odd heights to A 
particles, and the even heights to B particles. In the following, the term monolayer (ML) will 
denote one AB layer, i.e. the height difference between adjacent monolayers is 2. Inside the 
bulk of the crystal, there is an attractive interaction between the particles and their nearest 
neighbours which is isotropic in directions parallel to the surface. While there is no difference 
in the interaction of B particles between the bulk and the surface, A particles on the surface 
interact with the anisotropic interactions of the simplified version of our lattice gas model of 
the CdTe(OOl) surface introduced in section |2.2.4 . The other parameters of the model are 



defined in Figure 3.1. Thus, the Hamiltonian of the crystal becomes 



H = SABnAB + £B,bnB,b + £A,bnA,b + £hnA,h + £dnA,d + £xnA,x (3.1) 

Here, uab is the number of A-B bonds between the layers, n^^b is the number of nearest 
neighbour pairs of B atoms, nA,b the number of A nearest neighbour pairs inside the bulk, 
nA,h the number of surface A atoms next to a higher column and nA,d i''T'A,x) the number 
of bonds of surface A atoms to diagonal neighbours (nearest neighbours in x-direction) at 
the same height. The occupation of neighbouring sites in the y-direction with 74-atoms is 
forbidden, independent of the height difference between both atoms. Additionally, atoms 
must not stay in the wrong sublattice. 

In addition to diffusion to nearest neighbour sites we permit diffusion to diagonal neighbour 
sites. This process is important to obtain a reasonable mobility of A atoms on a i? terminated 
surface. An example is shown in figure p.2| a. In the absence of diagonal neighbour diffusion, 
due to the repulsive interaction in y-direction an A atom at the edge of an ^-cluster cannot 
diffuse along the cluster edge. Consequently, an equilibration of the cluster configuration is 
possible only via the detachment of atoms from the cluster. However, this process has a large 
activation energy since all bonds in the horizontal direction must be broken. 

Additionally, we permit the diffusion of an AB pair with the B atom on top of the A 
atom, if the diffusion of the B atom alone would end in the wrong sublattice. This process is 
required to preserve the ergodicity of the system. Consider the surface configuration shown 
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Figure 3.2: Motivation for diffusion to diagonal neighbour sites and diffusion of AB pairs. A 
atoms are shown as dark cubes, B particles in light grey. The x-axis points from left to right, 
the y-axis from front to back. Panel (a): Diffusion of A atoms on a 5 terminated surface. 
The diffusion of an A atom (arrow) along the edge of a cluster of other A atoms is severely 
restricted without diffusion to diagonal neighbour sites. Panel (b): Coalescence of islands is 
impossible without AB diffusion. 



in figure 3.2b. Since atoms must not stay in the wrong sublattice the B atoms on top of 
the islands cannot diffuse across the gap between the islands. Consequently, the energetically 
favourable coalescence of the islands is suppressed unless AB pairs are allowed to diffuse 
together. However, it is not necessary to permit also the diffusion of AB pairs with the A 
atom on top since the vacancy structure of A terminated islands is more flexible than a flat 
B terminated layer. We have performed some trial simulations where the diffusion of both 
arrangements of AB pairs was permitted. The results are identical to those of the standard 
model. 

Diffusion is an Arrhenius activated process (equation L2). For simplicity, we assume an 
equal activation energy Bq for all diffusion processes where the energy of the system does not 
increase. Then, due to the detailed balance condition B.3 which will be discussed in appendix 
P, the barrier height for the diffusion of a particle to a site which is energetically unfavourable 
compared to the initial site is necessarily Bq + AH. Here, AH is the energy difference between 
the initial and the final state of the system. In the literature, a dynamics with these rates 
is denoted as "Kawasaki dynamics" |Q. The desorption of a particle requires an additional 
activation energy e^. This dynamics does not depend on the parameter eab, so it needs not to 
be specified, and eB,b and eA,b enter only via the sum e' := eB,b + £A,b- The constant prefactor 
1/to := 1^ exp{—BQ/kT) sets the timescale of the model and needs not to be specified, if we 
measure time in units of to ■ Particles from an external source arrive at each lattice site with 
equal probability. If the adsorption of the particle at this site would lead to a forbidden state, 
it is reflected. 

The simulations are performed using continuous time Monte Carlo techniques. A detailed 
discussion of the algorithm can be found in appendix 

For simplicity, we measure energy in units of \ed\ and set Boltzmann's constant to unity. 
In the following, we use the parameter set e' = Ed = —I, £h = 0, £x = —1.97, = 1.5. With 
these parameters, the difference in the surface energies between perfectly c(2 x 2) and (2 x 1) 
reconstructed surfaces is small compared to the surface energies themselves. This is consistent 
with our finding in chapter |2| that a small AE is an essential feature of the c(2 x 2)cd-(2 x l)cd 
reordering. We have verified that the behaviour of our model remains qualitatively the same 
for a wide range of parameters, as long as the basic property of a small energy difference 
between the A terminated reconstructions is preserved. 
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3.2 Characterization of the surface 



The surface width W := {{hij — (/ijj))^) is a measure for the roughness of the surface. To 
characterize the arrangement of A atoms quantitatively, we investigate correlations and an 



order parameter similar to the quantities which were defined in section |2.2.2| . To this end, we 
consider an array c of integers Cij := hij mod2. Cij = 1 if the column of atoms at is 
terminated by an A atom and zero otherwise. These variables can be used to calculate the A 
coverage pA ■= N'"^ Z^f^=i Cij- and the normalized correlations 

1 ^ 

^ij (ci+i,j+i + Ci+i,i-i) (3.2) 



K'X 



2paN^ . . , 



— ^ ^ QjQ+i,,. (3.3) 



These definitions are similar to those of Ci,, and Cq,, in equations 2.2 and |2.3[ We have 



normalized the correlations with the A coverage such that K'^ and are equal to the 
fraction of A atoms which are incorporated in locally c(2 x 2) and (2 x 1) reconstructed 
regions of the surface, respectively. The long range order of the A atoms is measured by the 



(2x2) . . (2x2) 

order parameter the definition of which is identical to that of in equation |2^ 



As we will show, there is no long range order in the strict sense in our model such that 

(2x2) 

' is zero even at low temperature. However, there might be correlations between A 
atoms at intermediate distance. These are characterized by means of the autocorrelation 
function 

C{x) := {djCi+xj+y) - Pa (3-4) 

of c. If the positions of A atoms at a distance x = (x, y) are uncorrelated, we have C{x) = 0. 
Ornstein-Zernike theory ||T^ states that in thermal equilibrium the asymptotical behaviour of 
|C(x)| for large \x\ is 

\C{x)\ ~ Coo + A|x|-('^-i)/2 exp(-|x|/0 (3.5) 

where ^ is the correlation length and d is the dimension of the system. A nonzero value of 
Coo indicates long range order. This equation is not valid at the locus of a phase transition 
where the correlation length diverges. Note, that in a system with anisotropic interactions the 
constants Coo, ^ and ^ may have a directional dependence. If neighbouring sites are correlated, 
C{x) is always positive. If the values of Cij on neighbouring sites are anticorrelated, C{x) 
oscillates between positive and negative values. In our model, this is necessarily the case if x 
is parallel to the y-direction since there is a hard repulsion between A atoms on neighbouring 
sites in this direction. However, if x points in the x-direction, the configuration of neighbouring 
sites is anticorrelated only if the A atoms arrange preferentially in a c(2 x 2) pattern. On a 
(2 X 1) reconstructed surface neighbouring sites in the x-direction are correlated. If ~ K'^, 
these effects compensate each other. Consequently, on such a surface we expect a very small 
correlation length in the x-direction. 

3.3 Layer- by- layer sublimation 



Figure p.3| shows the evolution of an initially B terminated, flat surface under vacuum at 
a temperature T = 0.5. At the beginning of the simulation, pA increases from zero to an 
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Figure 3.3: Sublimation starting from a 
flat B terminated surface at T = 0.5. The 
system size is = 128. Solid line: W, 
dashed: KJ^, dotted: Kj^. For clarity of 
plotting, the last two curves have been 
smoothed; W shows the natural fluctua- 
tions. Each oscillation period corresponds 
to the desorption of one monolayer. 



asymptotic value close to 0.5. Thus, the sublimating surface is A terminated. 

We find, that this increase can be fitted with an exponential relaxation, pA_ = P^i^ ~ 
exp(— t/r)). Here, r is the time constant of the decay of the B terminated surface. The tem- 
perature dependence of the corresponding desorption probability p (x 1/t follows an Arrhenius 
law, p oc ex.p{—Esuh/T), with an activation energy E'sub = Bq + 3.0 ± 0.1 

After this onset, A and B evaporate stoichiometrically, and sublimation proceeds in layer- 



by-layer mode (see section 1.7). This is reflected by the oscillations of the surface width W 
(solid line in figure |3^ ) . Whenever the surface is atomically flat, W is minimal. The maxima 
of W correspond to a configuration where 50% of the surface are covered with islands of one 
monolayer thickness. The dependence of the sublimation rate on T follows an Arrhenius law 
with an activation energy Bq + 6.0 it 0.2 which is greater than that found for the decay of the 
B terminated surface. However, in our model the microscopic energy barriers for desorption 
of A and B are equal. Consequently, the difference in the macroscopic activation energies is 
solely a consequence of the stabilizing effect of the reconstruction of the A terminated surface. 

Surprisingly, K"^ and oscillate during layer-by-layer sublimation. Each time a com- 
plete AB layer has desorbed and the surface is atomically flat (minima of W) the c(2 x 2) 
reconstructed fraction of the surface is maximal. On the contrary, a rough surface with a large 
number of islands (maxima of W) seems to prefer the (2x1) reconstruction. This can be 
understood from the fact that the attractive lattice gas interactions are present only between 
particles in the same layer. Thus, the island edges impose open boundary conditions to the 
lattice gas of A atoms on the island. In contrast to a c(2 x 2) reconstructed domain, a (2 x 1) 
terminated island can reduce the energy of its boundary by elongating in x-direction. Since 
the ground state energies of both structures are nearly degenerate, the formation of (2 x 1) 
may reduce the surface free energy. This picture is confirmed by the fact that islands on 
sublimating surfaces are indeed elongated (see figure ^^b). 

3.4 Step-flow sublimation 

Clearly, an oscillatory behaviour cannot be understood within the framework of equilibrium 
thermodynamics. However, pure layer-by-layer sublimation is rarely ever observed in exper- 



iments. In particular, investigations of CdTe [73| have shown, that step-flow sublimation is 
the dominant mechanism. 

We determine the energy of steps in our model. This is done by calculating the energy 
which is needed to shift one half of a large crystal by the thickness of one monolayer in the 
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Figure 3.4: Evolution of the surface of a 64 x 64 system in layer-by-layer sublimation. A 
atoms are shown as filled circles, B atoms as open circles. Greater symbols denote greater 
surface heights. Panel (a) shows a surface snapshot at t = 7000 -to where the surface roughness 
is minimal, panel (b) shows the surface at t = 19000 • to where the roughness is maximal. 



z-direction. In this process, two steps are created: one at the upper side of the crystal and 
a second one at the lower side. Consequently, the step energy is one half of the total energy 
costs. Normalizing this with the length of the step, we obtain the step energy per unit length. 
On perfectly c(2 x 2) reconstructed surfaces, we find 

7diag = -^{ex + e') (3.6) 

for diagonal steps and 

7par = -£d - 2^ + (3.7) 

for steps oriented in x- or y-direction. With our standard parameter set, we have 7diag = 
a/2 < 7par = 3/2. Therefore, at T = steps oriented at an angle of 45° to the coordinate 
axes are favoured energetically. We expect, that this holds also at finite temperature as long 
as the surface is c(2 x 2) reconstructed. 

Diagonal steps in the simple cubic lattice of our model correspond to steps oriented in 
the [100] or [010] direction on the zinc-blende lattice. Indeed, investigations of CdTe(OOl) 



by means of STM microscopy [g^, 68] have shown, that this is the preferential orientation of 
steps. 

In the following, we present the results of simulations of vicinal surfaces with diagonal 
steps. We find that the oscillations in the correlations disappear for terrace widths smaller than 
~ 50 lattice constants, where sublimation proceeds in step-flow mode. Then, the correlations 
become stationary apart from statistical fluctuations. 

Figure 3^ shows the T-dependence of pA, K'^, K\ and the sublimation rate rgub in this 
stationary state at a step distance of 32\/2 lattice constants. In the investigated temperature 
range pA decreases only slightly from pA = 0.498 at T = 0.4 to pA = 0.41 at T = 0.8. Clearly, 
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Figure 3.5: Vicinal surfaces sublimating in step flow mode. The images show sections of 
64 X 64 lattice constants of surfaces at the end of the simulation runs shown in figure 3.6. 
Panel (a): T = 0.4 (below the reordering), Panel (b): T = 0.6 (above the reordering). 
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Figure 3.6: 

tions K^A 



^-coverage pA (+), correla- 
K'j^ ) and sublimation rate 
rsub (Q) in the stationary state of the sub- 
limation of vicinal surfaces under vacuum. 
The system size is = 128. The steps are 
oriented in the diagonal direction, the ter- 
race width is 32-v/2 lattice constants. The 
data points show time averages of the inves- 
tigated quantities. 3 • lO^to time has been 
simulated both for relaxation and for mea- 
surement. 



at r = 0.55 there is a transition from a c(2 x 2) configuration at low temperature to a high 
temperature regime where the (2 x 1) ordering dominates. This is reminiscent of the behaviour 
of the correlations in the vicinity of the order-disorder phase transition of the planar lattice gas 
(compare with figure |2.4| ). Figure |3.5| shows typical configurations of the surface at T = 0.4 
which is below the reordering (figure I 



i) and at T = 0.6 above the reordering (figure 3.5b). 
Comparing the results of our simulations of the sublimating surface with the anisotropic 



lattice gas introduced in section 2.2.4 we can check in how far the properties of the A atoms on 
the surface can be understood within an equilibrium theory of a flat surface. For simplicity, 
in the following we denote the planar lattice gas as "2D model" and the solid-on-solid model 
as "3D model". We investigate the 2D model at the value of pA which is measured on 
the surface of the 3D model sublimating in step-flow mode. We perform canonical Monte 
Carlo simulations where the number of A atoms is fixed. Additionally, we consider results of 
transfer matrix calculations where the chemical potential p is adjusted to obtain the desired 
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Figure 3.7: Panel (a): Correlations in the planar lattice gas at the value of pA which is 
found on a sublimating surface. The dashed (dotted) line shows (-f^l) obtained from a 
transfer matrix calculation at a strip width L = 16. Additionally, we show the results of 
canonical Monte Carlo simulations (x: K^, ^: Kj^. The system size is = 128. 10^ • N"^ 
Monte Carlo steps have been performed both for equilibration and for measurement. Panel 
(b): Correlation length in the y-direction measured on the sublimating surface (x) and in an 
equilibrium Monte Carlo simulation of the planar lattice gas (□). The solid line shows results 
of a transfer matrix calculation. These data have been measuered in the simulations shown 
in figure 3.6 and in panel (a). 



A coverage. Of course, the values of the parameters = —1.97 and = —1 are identical 
to those used in the 3D model. Figure 3^ shows the equilibrium values of and Kj^ 
in the 2D model. The simulation results agree well with the results of the transfer matrix 
calculation. Qualitatively, the behaviour of the 2D model resembles that of the 3D model. 
There is a crossover between a dominant c(2 x 2) arrangement at low temperature and a high 
temperature regime where the A atoms preferentially arrange in rows. However, there is no 
quantitative aggreement. In the 2D model, the reordering occurs at a higher temperature. 
Additionally, at all investigated temperatures the c(2 x 2) reconstructed fraction of the surface 
is greater than in the 3D model. 

We determine the temperature at which the 2D model leaves the long range ordered 
phase (I) and enters the disordered phase (II). Monte Carlo simulations and transfer matrix 
calculations yield Tc ~ 0.45 for the phase transition of the infinite system. However, in our 
simulations of the 3D model we find that = apart from statistical fluctuations even at 

T = 0.35. Thus, there is no long range order of the A atoms on the surface of the sublimating 
crystal. 

An inspection of surface snapshots of the 3D model (figure |3.5| a) at low temperature 
shows, that there are several c(2 x 2) reconstructed domains. In each domain, the A atoms 

f 2 X 2) 

occupy one of the two sublattices which correspond to positive and negative values of M)^ , 

f2x2) 

respectively. Domains of both orientations cover equal areas such that on average MX ' = 0. 
The domains have an elongated shape where the domain size in the y-direction is much greater 
than that in the x-direction. Frequently one encounters domains which run through the whole 
system. Thus, there are long domain walls parallel to the y-direction which consist of a row 
of pairs of A atoms on neighbouring sites in the x-direction. Of course, these pairs contribute 
to Kj^. Consequently, the formation of domain walls on sublimating surfaces explains both 
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the absence of long range order and the greater values of Kj^ compared to the 2D model. 

We investigate correlation lengths in the 2D and in the 3D model. This is done by 
fitting an Ornstein-Zernike law (equation |3.5| ) to C{x) which is measured in the Monte Carlo 
simulations. 

To check whether this method yields correct results we compare the correlation length in 
the y-direction in the 2D model with a transfer matrix calculation where is determined 



from the eigenvalues according to equation A. 16 . We find, that C((0, y)~^) follows an Ornstein- 



Zernike law for distances greater than 15 lattice constants. The results are shown in figure 



3.7|b. At high temperature, Monte Carlo data and results of the transfer matrix calculation 
agree well. Clearly, in a finite system with periodic boundary conditions, the correlation 
length cannot exceed one half of the system size. Therefore, at low temperature where the 
transfer matrix calculations yield very large ^y, the Monte Carlo results saturate at values 
close to 64. Consequently, our fitting procedure is reasonable for correlation lengths which 
are small compared to the system size. 

We find, that in the 3D model C((0, y)^) agrees well with an Ornstein-Zernike law although 
this is strictly valid only in thermal equilibrium. Interestingly, we find values of S^y which agree 
well with the equilibrium results of the 2D model (figure |3.7[ b) . 

Since C{{x,0)~^) depends strongly on which of and K'^ is greater, we expect a different 
behaviour of correlations in the x-direction in the 2D and in the 3D model. Indeed, this is 
the case. In general, correlation lengths in the x-direction are small. Consequently, on 
large lengthscales where a fit to an Ornstein-Zernike law would make sense, C{{x,0)~^) is 
indistinguishable from its asymptotical value. Therefore, we present only qualitative results. 
At all temperatures, in the 3D model is significantly greater than in the 2D model. Probably 
this is a consequence of the rows of A atoms parallel to the x-direction which are characteristic 
of the (2 X 1) arrangement. 

3.5 The surface under an external particle flux 

To get insight into the behaviour of a material in an MBE environment, it is important to 
understand how the structure of the surface changes if it is exposed to a particle beam. Ad- 
ditionally, these investigations give us insight into the behaviour of the surface at A coverages 
which differ from the value found under vacuum conditions. 

We simulate vicinal surfaces which are exposed to a flux of pure A or B. Panels (a) and 



(c) of figure 3^ show the quantities pA, and K'^ as functions of the particle flux F at a 
temperature T = 0.45 which is below the c(2 x 2)-(2 x 1) reordering under vacuum, and at 
T = 0.57 which is slightly above the transition. 

An A flux increases pA- This leads to an increase of K"^ such that under a large A flux the 
A atoms arrange preferentially in a c(2 x 2) pattern even at temperatures above the reordering 
under vacuum (figure |3.8| c). Applying B fluxes, we can regulate pA to values smaller than 
those close to 0.5 observed under vacuum. We find, that a decrease of pA makes the remaining 
A atoms arrange preferentially in the rows of the (2 x 1) reconstruction. 

In figure |3.8| b,d we show the correlations Kj^ and K"^ which are found in the 2D model 
at the value of pA which is measured in the 3D model. A comparison of these results with 
the behaviour of the 3D model shows, that the qualitative similarity and the quantitative 
differences between the 2D model and the 3D model which have been observed under vacuum 
are preserved in a wide range of A coverages. There is no long range order of the A atoms 



48 



(a) 



(b) 



0.75 



X * >K H 



0.75 



0.5 



0.5 



0.25 



0.25 



0.1 0.05 
<-B 



F-tn 



0.05 0.1 



0.1 



0.05 







0.05 



0.1 



A-> 



(c) 



0.75 



0.5 



0.25 y. 



(d) 



X X X >; 



X 

< 



0.75 



0.5 



0.25 



0.1 0.05 



F-tn 



0.05 0.1 



0.1 0.05 
^B 



F-tn 



0.05 0.1 



Figure 3.8: A coverage pA (+) and correlations Kj^ (x), K'^ (>|^) in the stationary state of 
vicinal surfaces under an external particle flux and equilibrium values in the planar lattice 
gas at the A coverage measured on the crystal surface. Panel (a): Behaviour of the 3D model 
at T = 0.45. The steps are oriented in the diagonal direction, the terrace width is 32^/2 
lattice constants. Panel (b): Results of the 2D model at the value of pA measured on the 
sublimating surface at the flux indicated at the rr-axis. The symbols denote results of Monte 
Carlo simulations, the lines show values determined by means of a transfer matrix calculation 
at a strip width L = 16. Panels (c) and (d): 3D model (c) and 2D model (d) at T = 0.57. All 
other parameters are identical to those used in (a) and (b). Note that an ^-flux applied to the 
surface of the 3D model leads to a re-entry into the c(2 x 2) reconstructed configuration at a 
temperature above the reordering under vacuum. In all Monte Carlo simulations shown, the 
system size is 128 x 128 lattice constants. In the equilibrium simulations, 10^- A^^ Monte Carlo 
steps have been performed both for equilibration and for measurement. In the majority of the 
simulations of the sublimating surface, time averages over 3 • lO^to; following a time interval 
of the same length for relaxation have been performed. Since the relaxation of surfaces under 
an A flux at T = 0.45 is very slow, a longer time interval of lO^to has been used in these 
simulation runs. 
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even at high A flux and low temperature. Instead, there are domains which are elongated in 
the y-direction. This yields ' = and values of which are greater than those found 

in the 2d model. At low A coverages we find that in the 3D model both and Kj^ are 
slightly smaller than in the 2D model. An investigation of C{x) at surfaces under a particle 
flux shows, that correlation lengths in the y-direction in the 2D model and in the 3D model 
agree quantitatively, while there are significant differences in 



3.6 Conclusions 

In summary, our simulations reproduce the experimental results for the reconstruction of 
Cd terminated CdTe(OOl) surfaces qualitatively. Under vacuum, there is a crossover from a 
c(2 X 2) reconstruction at low temperature to a high temperature phase with a dominant (2 x 1) 
arrangement. An external flux of A atoms stabilizes the c(2 x 2) arrangement at temperatures 
above the reordering in vacuum. We find that correlation lengths in the y-direction are much 
greater than those in the x-direction. This agrees with experimental results reported in 

We have shown that the reconstruction order of A atoms on a vicinal surface sublimating 
in step flow mode does not change qualitatively compared to a planar lattice gas in thermal 
equilibrium. In particular, the behaviour of the short range correlations K^^ and Kj^ in the 
3D model is quite similar to that found in the 2D model. This shows, that the investigation 
of comparatively simple models like those introduced in chapter ^ is an appropriate tool to 
get insight into the behaviour of reconstructed surfaces under MBE conditions. 

The nonequilibrium conditions of sublimation induce characteristic deviations of the ar- 
rangement of A atoms on the surface from the equilibrium conflguration. Perhaps the most 
important one is the absence of long range order at low T and high pA- This is reminiscent 



of the behaviour of the growth model of Chin and den Nijs [P0 [. 

Our findings can be understood within the following picture of sublimation: There is 
a permanent creation and annihilation of reconstructed surface. As the terminating layer 
of the crystal evaporates, the layer below is uncovered and becomes the new surface. In 
this process, a given atom stays on the surface for a typical time Trd = l/Tgub before it is 
desorbed. During this time, the arrangement of the atoms evolves towards the equilibrium 
configuration. However, it is insufficient for a complete equilibration. Therefore, there are 
characteristic defects in the reconstruction order of A atoms. 

At low temperature and high A coverage, there is a division of the surface in domains. 
Most of the domain walls run in the y-direction. This can be understood from the fact that the 
formation of such a domain wall requires only a small energy of /S.E = \2e(i — ea;|/2 = 0.015 
per unit length. Thus, there is only a comparatively small driving force which makes the 
domains coarse. Similarly, the small energy difference between a c(2 x 2) order and a (2 x 1) 
arrangement might be responsible for the high values of at high temperature. At small 
PA, both Kj!^ and K'^ are smaller than in thermal equilibrium. This indicates, that in the 
3D model there is a higher fraction of single A atoms than in the 2D model. This may be 
understood from the fact that Trei is not sufficient for some diffusing adatoms to meet a binding 
partner. 

On the other hand, the good agreement between the correlation lengths in the y-direction 
measured in the 2D model and in the 3D model shows, that important aspects of the statistical 
mechanics of the lattice gas in thermal equilibrium are preserved in step fiow sublimation. 
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Chapter 4 

Simulation of CdTe 



In this chapter, we extend the model of a binary compound with surface reconstructions pre- 
sented in chapter ^ to make it a more reahstic model of CdTe. We simulate the zinc-blende 
lattice of II- VI semiconductor crystals instead of a simple cubic lattice. To be meaningful, 
a fully three-dimensional modelling of the CdTe crystal requires a detailed understanding of 
the interactions of the atoms on all possible surface orientations and the correct treatment 
of vacancies and non-stoichiometric crystals. Such a model is beyond the scope of this in- 
vestigation which aims at an understanding of the behaviour of reconstructed (001) surfaces. 
However, as we will show, a physically relevant model is not possible within the limits of the 
solid-on-solid condition. 

To find a reasonable compromise between the complexity of a fully three-dimensional 
model and the limits of a strict solid-on-solid model, we permit the occupation of Cd sites 
by Te atoms and the binding of Te atoms to single Cd atoms on the surface, while the 
incorporation of point defects into the bulk of the crystal is forbidden. Then we can treat 
those atoms which stay in the correct sublattice and have at least two chemical bonds to 
atoms in the layer below in a solid-on-solid manner, while the weakly bound Te atoms are 
considered separately. 

First, we study the behaviour of this model in vacuum and under a flux of pure Cd or 
pure Te. In particular, we show that a reasonable choice of model parameters yields semi- 
quantitative agreement with experiments. Then, we perform simulations of atomic layer 
epitaxy and compare our results with those of experimental investigations of CdTe by means 
of this growth technique. 



4.1 The zinc-blende lattice 

In this section, we present a solid-on-solid representation of the zinc-blende lattice. A sketch 



of its unit cell is shown in figure |2.1| . We use cartesian coordinates where the x-axis points in 
the [110] direction, the y-axis points in the [110] direction and the z-axis points in the [001] 
direction which is perpendicular to the (001) surface. The origin of the coordinate system is 
at the centre of a Te atom. If the zinc-blende lattice is projected onto the x-y plane in the 
z-direction, the image positions of all atoms lie on a square lattice (figure [4.1| ). Measuring 
the lattice constant in appropriate units, the images of Te atoms are at postitions (2i, 2j) and 
(2z -|- l,2j -|- 1) in the x-y plane where i,j E Z. The images of Cd atoms are at positions 
(2i,2j -|- 1) and {2i + l,2j) with integer i, j. This square lattice is used as the basis of our 
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Figure 4.1: The zinc-blende lattice lattice seen from 
the [001] direction. The figure shows a Te terminated 
surface with a step in the x-direction. Cd atoms are 
shown as filled circles, Te atoms as open circles. The 
size of the atoms corresponds to their z-coordinate: the 
higher the atom, the larger the circle. 

solid-on-solid representation of the zinc-blende (001) surface, which is described by a two- 
dimensional 2L X 2N array h of integers, hx^y is the maximal z-coordinate of the column of 
atoms whose images are at {x,y), measured in units of the spacing between a Cd and a Te 
layer. In the following, for simplicity we will loosely denote any atom whose image is at (x, y) 
as "atom at {x,y)" and h^^y as "height of the atom at {x,y)" . 

An atom in the zinc-blende lattice may have at most 4 binding partners. Since the solid- 
on-solid condition forbids overhangs and vacancies, each atom must have at least two bonds 
to atoms of the opposite species in the layer below. Only atoms with exactly two bonds may 
be removed from the crystal without violating the solid-on-solid condition: An atom with 4 
bonds is inside the bulk of the crystal. If we remove it, we create a vacancy. If we remove an 
atom with three bonds the atom in the layer above which is bound to it becomes an overhang. 
Atoms with only one bond form an overhang and do not exist in a solid-on-solid model. In the 
following, we will denote atoms which may be removed from the surface without violating the 
solid-on-solid condition as mobile. Similarly, an atom can be deposited only at a site where it 
has exactly two bonds. The deposition of an atom with only one bond creates an overhang. 
Empty sites where the deposition of an atom creates three or four bonds can exist only in 
the presence of overhangs or vacancies, which is not the case in a solid-on-solid model. In 
consequence, a diffusion process where an atom is removed from one site and deposited at 
another cannot change the number of chemical bonds in the system. 

A Cd atom at (x, y) is bound to Te atoms at (x, y =t 1) in the layer below and can be bound 
to Te atoms at (x it 1, y) in the layer above. The atom is mobile, if there are no bonds to Te 
atoms in the layer above, which is the case if hx,y > hx±i^y The neighbouring Cd atoms of a 
Te atom at {x',y') are at {x' ± l,y') in the layer below and at {x',y' ± 1) in the layer above. 
Thus, the condition for mobility of a Te atom is hx'^y' > hx\y'±i- 

An atom can be deposited at a site if the solid-on-solid condition is fulfilled and the atomic 
species matches the character of the site. At a Cd site at {x,y), the solid-on-solid condition 
demands that hx,y-i = hx,y+i and hx,y < hx^y±i. A Te site at {x',y') can be occupied if 
hx'-i,y' = hx'+i,y' and /ia;'y < /ix'±i,y'- The deposition of an atom at {x,y) is performed by 
setting hx,y <— hx,y + 4. Similarly, setting hx^y <— — 4 removes the atom at {x,y). 

As discussed in section |2|, the positions of the atoms in one half-layer of the zinc-blende 
lattice lie on a square lattice. However, this is not identical to the square lattice structure 
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of the solid-on-solid representation of the crystal lattice, as it is in the simple cubic model 
presented in chapter ^. The distance between nearest neighbours in the square lattice of 
one half-layer (denoted as lateral neighbours in the following) is twice the distance between 
neighbouring images of atoms in the projection of the whole crystal onto the x-y plane. 

4.1.1 Desorption and adsorption 

Due to the solid-on-solid condition, only mobile atoms may desorb. We consider the desorption 
of single atoms from the surface, neglecting the fact that in reality Te probably desorbs as a 
binary molecule from the CdTe surface. Similarly, we consider only the adsorption of single 
atoms via an incorporation process (section 1.4.1| ): First, the species of the atom and the 



landing site are chosen randomly. Then, we search for sites where the depositon of this atom 
is possible within an incorporation radius around the landing site. The atom is deposited 
at the lowest site which has been found in that search. If two or more sites of equal height 
have been found, one of these sites is chosen randomly. If no appropriate site is found, the 
atom is reflected. 

4.1.2 Difrusion 

We consider an atom diffusing on a flat surface, for example a Cd atom on a Te terminated 
surface. For an atom at xq, the nearest sites that can be occupied by an atom of its species 
are the lateral neighbour sites at xq -\- Ax, where the hopping vector Ax is an element of the 
set {(0, 2), (0, —2), (2, 0), (—2, 0)}. Similar to the model discussed in section |3| we consider also 
diffusion to lateral diagonal neigbour sites, i.e. Ax G {(2, 2), (2, —2), (—2, 2), (—2, —2)}. 

However, a particle that can perform only these diffusion hops is not able to cross a step. 



We discuss this for the example of the Cd atom at site A in figure 4.1 diffusing downhill 
at a step in x-direction. It is not possible to deposit the atom at either of the sites B, C 
and D, since the condition h^^y-i = hx^y+i is not fulfilled there. Similar considerations hold 
for other step orientations, diffusion uphill and for the diffusion of Te atoms. Thus there 
is an infinite Schwoebel barrier in a model in which only these diffusion hops are consid- 
ered. Since we do not want to restrict ourselves in this special case, we introduce additional 



hopping vectors for diffusion. In figure [4.1| , the nearest sites to site A at the lower side of 
the step that can be occupied by a Cd atom are E and F, which yield Ax = (1,3) and 
Ax = (—1,3). Extending these considerations to different step orientations, we obtain the 
set {(3, 1), (3, -1), (-3, -1), (-3, 1), (1, 3), (-1, 3), (1, -3), (-1, -3)} of hopping vectors which 
are sufficient to permit the crossing of steps of one monolayer height difference oriented in 
arbitrary direction. To permit the crossing of steps higher than one monolayer, we would 
have to introduce additional hopping vectors. However, such high steps rarely occur in our 
simulations. Therefore, it is a reasonable approximation to neglect diffusion accross high steps 
in order to reduce the complexity of the model. In summary, there are 16 hopping vectors for 
the diffusion of an atom. However, due to the conditions that must be fulfilled to permit the 
deposition of the atom at the landing site, at most 8 diffusion hops are allowed for a single 
mobile atom. 
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4.2 The Hamiltonian of the solid-on-solid model 



The basic idea behind our model of the CdTe(OOl) surface is similar to that of the simple 
cubic crystal presented in chapter ^: besides the Cd-Te bond, there is an isotropic interaction 
between Cd atoms and between Te atoms on lateral nearest neighbour sites in the bulk. The 
interaction of Te atoms on the surface is identical to that of Te atoms in the bulk, i.e. we 
neglect the formation of Te dimers. The interaction of Cd atoms, however, depends strongly 
on whether they are mobile or incorporated into the crystal. Mobile Cd atoms interact with 
their lateral nearest neigbours and lateral diagonal neighbours with the interactions of the 
lattice gas model introduced in chapter |2[ 



4.2.1 The Cd-Te bond 

There is a chemical bond between Cd and Te atoms on nearest neighbour sites. We assume 
that the total energy of these bonds is proportional to their number with an energy Ec per 
bond. Since the strength of the interaction between atoms is a quickly decreasing function of 
their distance, we expect this energy to be the greatest part of the total binding energy of the 
crystal. To remove a mobile atom from the surface, two Cd-Te bonds have to be broken. This 
yields an energy barrier = —2ec for the desorption of an atom in addition to contributions 
from other, weaker interactions. 



As discussed in section 4.1, a diffusion process cannot change the total number of nearest 
neighbour bonds. Consequently, the strength of the Cd-Te bond does not influence the energy 
difference between the starting and the final state of the atom. In our model, the energy barrier 
a diffusing adatom must overcome is a separate parameter. It is smaller than the energy 
required to break up its bonds. Thus, the rates of diffusion processes are not a function of £c- 



4.2.2 Interactions with lateral neighbours 

Interactions with a longer range than nearest neighbour interactions are required to describe 
surface reconstructions and to ensure that the equilibrium surface at low temperature is flat. 

It is a reasonable approximation to consider only interactions between particles in the 
same layer. In a solid-on-solid model, the environment of an atom in the vertical direction is 
always the same. Therefore, interactions of longer range with atoms in this direction will be 
always present as long as the atom remains on the surface. In case of desorption they may be 
treated as an additional contribution to Ec- 

We introduce an attractive interaction between atoms which are incorporated in the crys- 
tal, i.e. atoms with more than two bonds. An atom interacts with its lateral nearest neighbours 
and the strength of the interaction is isotropic. However, we allow a different strength for the 
attraction between Cd atoms {£cd,b) and Te atoms (eTe,b)- 

We assume that the interaction between Te atoms is independent on whether they are 
incorporated in the crystal or mobile atoms on the surface and do not consider the dimerization 
of surface Te atoms. This simplification reduces the complexity of the model considerably. 
Mobile Cd atoms interact via the anisotropic interactions Ex, Ed and a hardcore repulsion 
between particles on lateral neighbour sites in the y-direction which were introduced in chapter 
^ Additionally, we introduce an interaction Ecd,h between mobile and incorporated Cd atoms 
which we assume to be isotropic. 
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Figure 4.2: Deposition of a Te atom on a c(2 x 2)(;;;d reconstructed surface via a Te* state. 
The Te* atom is shown as an open circle with a star inside. 

In summary, the Hamiltonian of the sohd-on-sohd model is 

Hsos = £cnc + £cd,bncd,b + £Te,bnTe,b + £:hncd,h + £xncd,x + £dncd,d (4.1) 

^ V ' 

Hlat 

Here, Uc is the number of chemical bonds beween Cd and Te atoms on nearest neighbour 
sites. ncd,b a-^d nTe,b a-re the numbers of pairs of incorporated Cd and Te atoms on lateral 
nearest neighbour sites. ncd,h is the number of mobile Cd atoms which are lateral nearest 
neighbours of incorporated Cd atoms. ncd,x 

and ncd,d are the numbers of pairs of mobile Cd 
atoms on lateral nearest neighbour sites in the x-direction and on lateral diagonal neighbour 
sites, respectively. In addition to these interactions we impose the constraint that the number 
^Cd,y of mobile Cd atoms on nearest neighbour sites in the y-direction must always be zero. 



4.3 Beyond the solid-on-solid condition 

In the present state of our model, there is a violation of the ergodicity condition similar 



to the one discussed in section 3.1. Consider the surface configuration shown in figure 4.3a. 
Energetically, it is favourable to shift the small island by one lattice constant into the negative 
y-direction to obtain the configuration shown in figure [4.3| e. However, it is not possible to 
construct a sequence of allowed single particle moves to do so. In chapter ^ we introduced 
the diff'usion of AB molecules with B atoms on top to solve an analogous problem. The 
counterpart of this process on the zinc-blende lattice would be the diffusion of a cluster of 



two Cd atoms and one Te atom. Then, the small island in figure [4.3| could diffuse as a whole. 
However, this is not a general solution to the ergodicity problem. It is not possible to remove 
such a cluster from the edge of a greater Te terminated island without violating the solid-on- 
solid condition. There is at least one additional Te atom in the top layer of the island which is 
bound to a Cd atom in the cluster. If we remove the cluster, this atom becomes an overhang. 

Additionally, in the present form of our model it is not possible to deposit an atom on a 
perfectly c(2 x 2)cd reconstructed surface. There are no sites which can be occupied by a Te 
atom without violating the solid-on-solid condition. The hardcore repulsion between mobile 
Cd atoms in the y-direction forbids a completion of the topmost Cd layer. 

There is a solution to both of these problems which has the advantage of being consistent 
with an experimental observation made on CdTe: under a strong flux of Tellurium one finds 
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Figure 4.3: Coalescence of islands. The Te atom on top of the small island stays temporarily 
in a Te* state. 



surfaces with a Te coverage of 3/2 and a (2 x 1) reconstruction (figure p.3[ ). This is explained 
by the formation of Te trimers which arrange in rows in the x-direction. A trimer is created if 
a Cd site on the surface is occupied by a Te atom. We consider this arrangement in our model. 
Additionally, we allow Te atoms which have only one bond to a Cd atom. In the following, we 
denote Te atoms in one of these configurations which both violate the solid-on-solid condition 
as Te* atoms. It is not necessary to give up a solid-on-solid model in favour of a fully three- 
dimensional model to consider Te* states if Te* atoms may exist on the surface of the crystal 
but the incorporation of point defects into the crystal is forbidden. We describe the crystal 
by two two-dimensional 2L x 2N arrays of integers, h and t. h represents those atoms which 



fulfil the solid-on-solid condition in the manner described in section 4.1. t represents the Te* 
atoms. A Te* atom at {x,y) corresponds to tx^y = 1, otherwise t^^y = 0. 
The rules for Te* atoms are the following: 

1. Te* atoms are allowed on sites with a mobile Cd atom (Te atoms with one bond to a Cd 
atom) and on sites where the deposition of a Cd atom would not violate the solid-on- 
solid condition (middle Te atom in a trimer). The coordinates of these sites correspond 
to Cd sites in the solid-on-solid model. Each site may be occupied by at most one Te* 
atom. 

2. Te* atoms on different sites do not interact. The binding energy e* of a Te* is indepen- 
dent on whether it is bound to two Te atoms or to one Cd atom. 

3. Two laterally neighbouring sites in the y-direction may be occupied by two Cd atoms, 
if at least one of the sites is occupied by a Te* atom. Thus, a Te* atom neutralizes the 
hardcore repulsion between laterally neighbouring Cd atoms in the y-direction. 
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4. Atoms on the solid-on-solid surface which are mobile may diffuse or desorb even if they 
are binding partners of Te* atoms. If a Te* atom is on a site which may not be occupied 
by a Te* atom any more after such a process (according to rule (1)), it is desorbed. Of 
course, the activation energy for the desorption of Te* atoms has to be considered in 
the calculation of the rate of the process. Te* atoms do not follow the motion of their 
binding partners. 

5. Both the desorption of Te* atoms and the deposition of Te atoms from an external 
particle beam in Te* states are possible. If a Te atom is landing on the surface, in 
the incorporation process Te* states within are considered in addition to states in 
the solid-on-solid lattice. Te* atoms may diffuse on the surface. We permit the same 
hopping vectors as for the diffusion of atoms on the solid-on-solid lattice. 

6. Te atoms may change between states in the solid-on-solid lattice and Te* states. This is 
a diffusion process where the Te atom hops between Te sites in the solid-on-solid lattice 
and Te* sites, which have the coordinates of Cd sites in the solid-on-solid model. We 
consider the set 

{(1,0), (-1,0), 

(1, 2), (-1, 2), (1, -2), (-1, -2), (2, 1), (2, -1), (-2, 1), (-2, -1), 
(3,0), (-3,0), (0,3), (0,-3)} 

of hopping vectors for this process. 

Thus, the Hamiltonian of our model is 

H = Hiat + ecn, + e*n\ (4.2) 

Here, Hiat is the contribution of the lateral interactions to the energy of the system. It has 



been defined in equation 4.1. ric is the number of Cd-Te bonds in the solid-on-solid lattice, 
n* is the number of Te* atoms. Additionally, we consider only surface configurations where 
the number n^^ ^ of pairs of mobile Cd atoms without a Te* atom on lateral neighbour sites 
in the y-direction is zero. 

Now we show examples of the role of Te* states in the deposition of a Te atom on a Cd 



terminated surface and the coalescence of islands. Figure |4.2| shows the deposition of a Te 
atom on a c(2 x 2)cd reconstructed surface. First, the atom is deposited in a Te* state in one 
of the gaps of the c(2 x 2)cd reconstruction (figure [4.2^ ). According to rule (3), one of the 
neighbouring Cd atoms can diffuse into the gap (figure [4.2|b). Following rule (4), now the Te* 
atom is bound to that Cd atom. This process creates a site in the solid-on-solid lattice which 
can be occupied by a Te atom. In the last step, the Te* atom diffuses to that site (figure 
f4.2| c) following rule (6). Figure shows the coalescence of islands. The Te atom on top of 
the small island stays temporarily in a Te* state. Then, the two Cd atoms can diffuse to the 
step edge. Finally, the Te atom is deposited on the site in the solid-on-solid lattice provided 
by the Cd atoms. 

The assumption of noninteracting Te* atoms (rule (2)) is a simplification which has been 
introduced to make an investigation of this model computationally feasible. Considering 
interactions between the quickly diffusing Te* atoms would lead to a considerable slowing 
down of the simulations. This simplification is justified if the density of Te* is so low that the 
probability to find a partner to interact with is small. We will show that interactions between 
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Figure 4.4: The potential energy surface of adatoms 



Te* atoms sliould be considered to describe the behaviour of a surface under a strong flux of 
pure Te at low temperature. 

A Te* atom which is bound to a Cd atom has three dangling bonds. As discussed in section 
1.1, these bonds are negatively charged. It is plausible to assume that this compensates the 
energetically unfavourable concentration of positively charged Cd dangling bonds on a pair 
of neighbouring Cd atoms in the y-direction. This might justify rule (3) which has been 
introduced to make the landing of Te atoms on a perfectly c(2 x 2)c(j reconstructed surface 
possible. Without rule (3), it would be impossible for a Cd atom to perform a diffusion process 
which provides a site in the solid-on-solid lattice that can be occupied by a Te atom. 

In the set of hopping vectors introduced in rule (6), we have considered diffusion to those 
sites which have a smaller distance to the starting site than the final sites of a diffusion process 
within the same sublattice. These hopping vectors permit the crossing of steps with a height 
of one monolayer. We have omitted the hopping vectors (0, 1) and (0, —1), since due to the 
topology of the zinc-blende lattice it is never possible to deposit the atom on the final site of 
such a diffusion process. 



4.4 The dynamics of adatoms 

In our model, we assume Arrhenius rates (equation ^]^) for diffusion and desorption processes. 
To calculate the activation energy we assume the energy of the transition state between the 
starting and the final configuration of the system to be equal to the greater of the energies 
of both configurations plus an additional barrier height. Since atoms on the solid-on-solid 
surface are bound stronger than Te* atoms, it is reasonable to assume that the barrier height 
for diffusion of atoms on the solid-on-solid surface is greater than that for diffusion of Te* 
atoms. In our model, we use the simplest choice for the energy barriers of all processes that 
meets this condition. There is an equal barrier height Bq for all diffusion hops of particles 
without lateral nearest neighbours and a barrier height B* for the diffusion of Te* atoms. 
The barrier heights for the other processes are chosen such that the barriers for all moves of 
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one particle which do not increase the energy of the system are equal and detailed balance is 
fulfilled for diffusion processes. 

This principle yields the following activation energies for diffusion on the solid-on-solid 
surface {Eg^g), diffusion of Te* atoms diffusion of a Te atom from a solid-on-solid state 

to a Te* state {Es_,^), diffusion of a Te* atom into a solid-on-solid state (E^^s), desorption 
of an atom on the solid-on-solid surface (Eg—,^) and desorption of a Te* atom 

Es^s = ma^{BQ,Bo + AH} (4.3) 

E^^^ = ma^{B* , B* + AH} (4.4) 

Es^* = maii{Bo,B* + AH} (4.5) 

E^^s = max{B* , Bo + AH} (4.6) 

Es^^ = ma^{Bo,AH} (4.7) 

E^^^, = max{B*,AH}. (4.8) 



Figure 4.4 shows a sketch of the energy barriers in our model. 

In the following, we will focus on the parameter set rj = 2, Ex = —1.95\ed\, Sh = 0, 
£cd,b = £Te,b = -0.8\ed\, Ec = -6.5|ed|, e* = -3|erf|, Bq = 9|ed| and B* = 2\ed\. Since there 
is an attractive interaction between mobile Cd atoms on lateral diagonal neighbour sites, 
£d is necessarily negative. This choice of parameters is guided by the following principles: 
Ex = — 1.95|ed| ensures that c(2 x 2)cd is energetically favourable compared to (2 x l)cd and 
the difference between the surface energies of both reconstructions is small compared to the 
total surface energy. The parameters £cd,b-, ^Te,b and £h have been chosen such that we obtain 
a c(2 X 2)cd reconstruction under vacuum at low temperature. Our model of a flat surface 
in thermal equilibrium (chapter ^) indicates that \£d\ is slightly smaller than 0.1 eV. On the 
other hand, we expect the activation energy 2£c which is required for the desorption of a 
single adatom to be on the same order of magnitude as the macroscopic activation energy 
of sublimation. Experimentally, values between 0.96 eV and 1.95 eV have been measured 



|21, 73, 75| depending on particle species and mode of sublimation. Thus, £c = — 6.5|erf| 



seems to be a good guess. An adsorption of Te on a c(2 x 2)cd reconstructed surface via a 
Te* state is possible only if Te* atoms are not desorbed too fast such that e* must not be too 
small. The barrier for diffusion of a Te atom on the solid-on-solid surface should be smaller 
than that for diffusion into a Te* state such that Bq < \ 2£c — e* | . Finally, the diffusion barrier 
for a Te* atom must be smaller than that for desorption of Te* which implies B* < |e*|. The 
prefactor of the Arrhenius rates has been set to = 10^^ /s which has turned out to be a good 
choice for CdTe [100|]. 

The value of \£d\ sets the energy scale of the model which needs not to be specified if 
temperature is expressed by means of the dimensionless parameter Q = kT/\£d\. If we set 
k = \£d\ = 1 like in the previous chapters, we have = T. In this chapter, we will determine 
\£d\ in physical units by identifying the temperature of the c(2 x 2)cd-(2 x l)cd transition with 
the experimental value. As we will show, this yields activation energies for sublimation which 
are on a reasonable order of magnitude. 

4.5 Characterization of the surface 

In principle, the methods applied for the characterization of the surface of this model follow 
the ideas introduced in chapters ^ and|3|. Note however, that there is an important difference 
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between the simple cubic lattice and the zinc-blende lattice: In our representation of the zinc- 
blende lattice, the distance between lateral nearest neighbours in one half-layer is two unit 
lengths. The positions of atoms in each CdTe monolayer are shifted by one unit length both 
in the x- and in the y-direction compared to the atoms in the layer below or above. Thus, 
the pattern of a surface reconstruction cannot be continued across a step as it is possible 
in the simple cubic lattice. Consequently, it is not possible to define the order parameter 



(2x2) 
Cd 



of the whole system in analogy to the planar lattice gas. Similarly, this effect hinders 



M, 

a calculation of correlation lengths with the method introduced in section 3.2. Therefore, we 
restrict ourselves to an investigation of the short-range correlations and K^^. 
To calculate these, we define a 2L x 2N array of integers c, where c,; 



1 if there is a 



mobile Cd atom at and zero otherwise. In terms of these variables, we have 



Cd 



^ 2L,2N 
2pcdNL ^ C,,,(q+2,,+2 + Q+2,,-2) 



K: 



1 



Cd 



where pcd = {LN) ^ Yl'i^=i the density of mobile Cd atoms on the surface. 



PCdNL 



2L.2N 

E 



(4.9) 
(4.10) 



4.6 Sublimation 

Qualitatively, the behaviour of a surface under vacuum is quite similar to that observed in 
our model of a simple cubic crystal with surface reconstructions presented in section ^. 

On a flat surface sublimating in layer-by-layer mode, the correlations and os- 
cillate (figure |4.5| ). Like in the simple cubic model, island boundaries induce a preferential 
arrangement of Cd atoms in rows such that there is a maximum of whenever the surface 
width W is large and a maximum of Kq^ whenever the surface is flat. 




Figure 4.5: Surface width l^(solid) and cor- 
relations K^^ (dashed), K^^ (dotted) during 
the sublimation of a flat surface at = 0.525 
which corresponds to T = 525 K in physical 
units. For clarity of plotting, the curves of 
the correlations have been smoothed; W shows 
the natural fluctuations. The system size is 
L = N = 6A. 



Figure ^^ a shows the correlations K^,^, K^,^ and the Cd coverage pcd as functions of 
temperature in the stationary state of step flow sublimation. The steps are parallel to the 
[100] direction and have a distance of 16 lattice constants. In the investigated temperature 
range, pcd is close to 0.5. It decreases slightly from pcd = 0.49 at O = 0.45 to pcd = 0.43 
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Figure 4.6: Panel (a): Stationary values of pcd {+), -f^cd (^) ^^"^ ^Cd (^) ™ ^*^P 
sublimation. Shown are time averages of these quantities which have been measured during 
the sublimation of 5 monolayers. Additionally, averages over 3 independent simulation runs 
have been performed. Panel (b): Equilibrium values of these quantities in the planar lattice 
gas in thermal equilibrium. The symbols show results of canonical Monte Carlo simulations 
at the values of pcd which are found on sublimating surfaces. 10^ • LN events have been 
performed both for equilibration and for measurement. The lines show results of a transfer 
matrix calculation at a strip width L = 16. The chemical potential /x has been adjusted 
to obtain the desired value of pcd- In both Monte Carlo simulations, the system size is 
L = N = 64. 



at = 0.7. At low temperature, the Cd atoms arrange preferentially in a c(2 x 2)cd order. 
At high temperature, the Cd atoms prefer to arrange in rows which are characteristic of 
the (2 X l)cd reconstruction. The crossover between both regimes where = is at 
G = 0.57. 

We compare the values of and which are found on the sublimating surface with 
those measured in the planar lattice gas introduced in section 2.2.4 at the same temperature 
and Cd coverage. We find, that the nonequilibrium conditions of sublimation enhance the 
tendency of Cd atoms to arrange in rows. At all temperatures, in step flow sublimation we 
find greater (smaller) values of i-^cd) than on a flat surface in thermal equilibrium. An 
inspection of surface snapshots shows, that at low temperature the c(2 x 2)cd reconstructed 
surface is divided into domains which are elongated in the y-direction. Thus, there is no 
long range order in the lattice gas of Cd atoms on one terrace. These observations parallel 
our findings in the simple cubic model. Consequently, the ideas about reconstruction order 
on sublimating surfaces developed in section ^ are valid also in the more complicated model 
presented in this chapter. 

We set the energy scale of our model by identifying G = 0.57 where Kq^^ = K^^ with the 
experimental value of w 570 K pl| , |75| for the temperature of the transition between c(2 x 2)^(1 
and (2 x l)cd- Then, we have T = G • 1000 K which corresponds to £d = — 0.0862 eV. 

Given the parameter set in physical units, we check its physical relevance by comparing 
activation energies of sublimation processes with the corresponding experimental values. A 
Te terminated surface decays in a temperature-dependent time interval r which corresponds 
to a sublimation rate rxe = 1/t. The measurement of sublimation rates in stoichiometrical 
evaporation is done straightforward by counting the number of desorbed particles. It is 
convenient to plot the logarithm of the sublimation rate versus 1/T. In this representation 
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which is denoted as Arrhenius plot, the data points he on a straight hne if the dependence of 
the sublimation rate on temperature is given by an Arrhenius law. 



Figure 4.7 shows that this is the case both for rxe,stcp ^ind for the rate rgtep at which CdTe 
evaporates stoichiometrically in step flow mode. The corresponding activation energies are 
-E'Te.stcp = 1.23 eV and -Egtcp = 1.43 eV. Experimentally, values of i^Tc.step = 0.96 eV [119] 



and -Bstcp = 1-54 eV |75] have been found. In layer-by-layer sublimation, the agreement of 
the data with an Arrhenius law is as good as in step flow sublimation. We find an activation 
energy .Biayer = 1.50 eV which is smaller than the experimental value 1.94 eV |73]. The 
activation energy for the decay of the Te-terminated surface in layer-by-layer sublimation is 
-^'Te.iaycr = 1.35 cV. Surprisingly the difference between Estep and -Biayer is quite small. From 
the quality of the fits of simulation data with an Arrhenius law, we estimate the accuracy 
of the values of the activation energies to about ±0.02 eV. For the attempt frequencies, we 
obtain values on the order of magnitude of 10^^/s. In general, the deviation of our values of 
the activation energies from the experimental results is less than 30%. This shows that our 
model parameters are at least in a physically reasonable region of the parameter space. 

At present, an optimization of model parameters to obtain better agreement with ex- 
periments does not make much sense. First, the experimental data are insufficient for a 
systematical fit of the nine parameters in our model and there are no values of diffusion barri- 
ers available from ab initio or at least semi-empirical calculations. Therefore, any parameter 
set will be partly based on estimates and physical reasoning. Second, our model still contains 
several important simplifications like a neglecting of Te dimerization and interactions between 
Te* atoms. Additionally, we have made a comparatively simple ansatz for the potential energy 
surface. Therefore, a perfect agreement with experiments is not expected. 
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Figure 4.7: Arrhenius plot of the sublima- 
tion rate in step-flow mode (+) and the rate 
at which the initial Te terminated surface de- 
cays (x). These data have been obtained in 
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the simulation runs shown in figure 4.6 a. The 
lines show fits with an Arrhenius law. 



4.7 Surfaces under an external particle flux 

Figure |4.8| a shows the Cd coverage and the correlations K^^, on a stepped surface which 
is exposed to a flux of pure Cd or pure Te at T = 500 K. Under vacuum, the Cd atoms 
arrange preferentially in a c(2 x 2)c(j pattern at this temperature. A Te flux decreases the Cd 
coverage of the surface. The remaining Cd atoms arrange preferentially in rows. Thus, the 
behaviour of our model at low pcd agrees qualitatively with that of the planar lattice gas in 
thermal equilibrium. Of course, this is what we expect in analogy to the simple cubic model 
introduced in section ^. 
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Figure 4.8: pcd (+)> -f^cd (^) ^Cd (^) stepped surfaces with a terrace width of 
16 lattice constants under an external particle flux. The system size is L = A'^ = 64. Panel 
(a): T = 500 K. The simulation shows time averages over 140s of simulated time after a 
relaxation phase of 140s. Panel (b): T = 600 K. The surface is exposed to a flux of pure 
Cd. Time averages have been performed over 7 s of simulated time after a relaxation phase 
of equal length. 



However, a Cd flux does not increase the tendency of the Cd atoms to arrange in a 
checkerboard pattern as it is the case for the A atoms in the simple cubic model. For Cd 
fluxes smaller than ~ 0.6 ML/s, the correlations Kq^ and Kq^ remain almost constant. At 
greater Cd flux, K^^^ increases and i^^^ decreases. A similar behaviour is found at T = 600 K 
where K^^^ > K^^ under vacuum (figure |4.8| b). For Cd fluxes smaller than 10 ML/s, K^,^ 
and remain close to their vacuum values. A greater Cd flux leads to an increase of Kq^. 

Surprisingly, there is no significant difference in the Cd coverage between a surface in 
vacuum and a surface which is exposed to an external Cd fiux. If we determine the equilibrium 
values of Kq^ and Kq^ in the two-dimensional lattice gas at the values of pcd which are 
measured on a surface under a Cd fiux, we obtain values which are almost identical to those 
found at the pcd of a surface in vacuum. This agrees qualitatively with the result that on 
a surface of a three-dimensional crystal which is exposed to a small Cd fiux and Kq^ 
are close to the values which are measured in vacuum. However, the increase of under a 
large Cd flux cannot be understood within the framework of a qualitative agreement of our 
model with the two-dimensional lattice gas. 

On the zinc-blende lattice, Cd atoms can be deposited only at sites where they are bound 
to two Te atoms in the layer below. This implies, that the maximum possible density of 
mobile Cd atoms on a stepped surface is smaller than 0.5. The exact value depends on the 
shape of the steps. Therefore, the fact that we measure pcd < 0.5 does not imply that there 
are many sites where an additional Cd atom can be deposited. If the density of such sites 
is small, an additional Cd flux cannot increase pcd further such that there is no reason for 
the mobile Cd atoms on the surface to arrange preferentially in a checkerboard pattern. This 
might explain why the behaviour of our model deviates from experiments which show that a 
Cd flux stabilizes the c(2 x 2)c(j reconstruction at high temperature. 

We expect a different behaviour if the dimerization of Te atoms in the terminating layer is 
considered. In a dimer, there is a chemical bond between two Te atoms which might also be 
the binding partners of a Cd atom in the layer above. Consequently, dimers are sites where 
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Figure 4.9: Sections of 16\/2 x 16^/2 lattice constants of surfaces which have been exposed to 
a flux of pure Cd for 280 s at T = 500 K. Panel (a): F = 0.6ML/s. Panel (b): F = lOML/s. 
Cd atoms are shown as filled circles, Te atoms as open circles. Greater circles denote greater 
values of hij. 



a Cd atom from a particle beam can be deposited. The binding energy of a dimer stabilizes 
the Te atoms with respect to desorption. Therefore, we expect that there are dimers on a 
surface sublimating under vacuum at high temperature if the binding energy of dimers is large 
enough. Such a surface is probably similar to the disordered phase of the planar lattice gas 
with Te dimerization introduced in chapter |2| (figure |2.7| e). Under a Cd flux, Cd atoms are 
incorporated into the dimers which leads to an increase of pcd- This makes the Cd atoms 
arrange preferentially in a c(2 x 2)cd pattern. 

The morphology of a surface under a small Cd flux is similar to that of a surface in 
vacuum. An example is shown in figure l^a. At large Cd flux, the morphology changes. In 
figure 4.Sb, we show a snapshot of a surface which is exposed to a Cd flux of 10 ML/s at 
T = 500 K. The terraces are not flat. Instead, there is a network of grooves with a depth of 
one CdTe monolayer which run across the surface. The edges of the grooves are (111), (HI), 
(111) and (111) facets. Between the grooves, there are ridges which consist of long rows of 
Cd atoms which are parallel to the x-direction. These rows yield the extremely large values 
of ~ 0.8 which are measured on such a surface. Consequently, the increase of K^^ on 
surfaces under a strong flux of Cd is induced by a facetting of the crystal surface. 

At low temperature, the increase of which encompanies the formation of grooves 
and ridges occurs at a smaller flux (0.6 ML/s at 500 K) than at high temperature (10 ML/s 
at 600 K). This indicates that under a small flux facetting is suppressed by the thermally 
activated diffusion of particles which smoothens the surface. A quantitative estimate confirms 
this consideration. We assume that the onset of facetting occurs at a fixed ratio between the 
Cd flux and the rate at which adatoms diffuse. In our model, the latter is given by an Arrhenius 
law with the energy barrier defined in equation 4.3 . Attractive interactions between particles 
in one layer are small compared to Bq. Therefore, Bq is a good approximation for the typical 
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energy barrier of a diffusing atom. This yields a ratio of 20 for the fluxes where facetting 
begins at T = 600 K and T = 500 K which agress well with the results of our simulations. 

The facetting of surfaces of II- VI compounds under an excess of one constituent has also 
been observed in experiments. In |129| ] it is reported that the (001) surface of ZnSe is facetting 
if it is exposed to a mixed particle beam with a large ratio of Zn to Se. The critical ratio 
between Zn and Se increases with temperature. 



4.8 Atomic layer epitaxy 

Atomic layer epitaxy (ALE) is used frequently in the fabrication of semiconductor heterostruc- 
tures and has been applied successfully to determine the stoichiometry of the surface recon- 
structions of CdTe and other II- VI semiconductors. Thus, it provides a scenario for the 
investigation of our model in a situation present in technical applications and permits a di- 
rect comparison of computer simulations with experiments. As explained in the introductory 



section 1.6, the basic idea of ALE is to obtain self-regulated growth by alternate deposition 
of pure Cd and Te. In the absence of reconstructions, one would expect the formation of 
complete half-layers of Cd and Te during deposition of the elements such that in each cycle 
one CdTe monolayer is grown. However, this is not necessarily the case in the presence of 
surface reconstructions where the number of atoms in the terminating layer of the crystal is 
different from that in a layer inside the bulk. 

In CdTe experimental investigations of the ALE growth rate Vg, i.e. the number of CdTe 
monolayers deposited per cycle have shown that it decreases in subsequent plateaus as tem- 



perature is increased |21, 27, 35 1. Within large temperature intervals, Vg is almost constant. 
At temperatures below ~ 510K, CdTe grows at a speed of IML/cycle. In the temperature 
interval 530 K < T < 640 K growth rates which are slightly lower than 0.5ML/cycle have 
been found. At T > TOOK, ALE growth is not possible at all. 




t(s) 



Figure 4.10: Thickness of the deposited layer 
(solid), K^^ (dashed) and K^^ (dotted) as 
functions of time in atomic layer epitaxy at 
T = 500 K. Pulse time and dead time are 
tp = Is and td = 0.1s. The particle flux dur- 
ing the pulses is 5 ML/s. The system size is 
L = N = 128. For clarity of plotting, the 
curves of K^, and K^^, have been smoothed. 



Since experiments [123| indicate that ALE growth proceeds in layer-by-layer mode, we 
focus on simulations of flat surfaces. Figure 4.10| a shows the evolution of pcd, and 
during ALE at T = 500 K. The initial configuration is a flat Te terminated surface. In the 
first phase, a Cd flux of 5 ML/s is applied for a pulse time = Is, followed by a dead time 
trf = 0.1s. In the second phase, Te is deposited with the same flux and identical parameters 
tp and td- This cycle is repeated permanently. Snapshots of the surface during ALE growth 
are shown in figure 4.11| . 
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Figure 4.11: Sections of 16\/2 x 16\/2 
lattice constants of the surface during the 
simulation run shown in figure 4.1C. Panel 
(a): t = 1.1s. Panel (b): t = 2.2s. Panel 
(c): t = 3.3 s. Panel (d): t = 4.4 s. Panel 
(e): t = 5.5 s. Cd atoms are shown as filled 
circles, Te atoms as open circles. Note the 
alternation of a flat surface with a surface 
covered with islands in subsequent cycles. 
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Figure 4.12: Temperature dependence of the 
growth rate in ALE on a flat surface. Vg has 
been determined from the increase of {hij) in 
cycles 2-4. Apart from T, all parameters are 



T(K) identical to the simulation shown in figure 4.1C. 



The first Cd pulse creates a Cd terminated surface with pcd ~ 1/2 (figure |4.11| a). Since T 
is below the c(2 x 2)c(j-(2 x l)cd transition, the dominant reconstruction is c(2 x 2)q^. At the 
onset of the Te phase, pcd decreases rapidly. This leads to an increase of K'^^ since at low 
Cd coverages it is energetically favourable for the remaining Cd atoms to arrange in rows. At 
the end of the cycle, approximately 50% of the surface are covered with Te terminated islands 
(figure 4.11b). The formation of a closed layer is impossible since only one half of the Cd atoms 
needed for a complete monolayer are present on the surface. The following Cd phase deposits 
half a monolayer of Cd again (figure 4.11| c). However, now the reconstruction is preferentially 
(2 X l)cd- This is due to the influence of the island edges on the reconstruction and analogous 
to the maximum of i^^'d observe in layer-by-layer sublimation after half a monolayer has 
desorbed. During the following Te phase, Cd atoms from the terminating layer of the islands 
and newly deposited Te atoms diffuse into the gaps between the islands which leads to the 
formation of an almost complete Te terminated surface (figure [4.11| d). In the next Cd phase, 
a c(2 X 2)cd reconstructed flat surface is created and the whole process repeats (figure 4.11| e). 
Thus, at the end of the odd cycles, the surface is rough, i.e. covered with islands, whereas at 
the end of the even cycles the surface is atomically flat. In our simulations the alternation of 
a rough and a flat surface can be traced during 3-4 cycles. There are two effects which keep 
it from repeating infinitely. First, we obtain a growth rate of 0.45 ML/cycle which is smaller 
than the ideal value of 0.5 ML/cycle. Second, the diffusion of particles from the islands into 
the gaps is not complete. Both effects hinder a perfect closure of a monolayer at the end of 
the even cycles which leads to a damping of the characteristic behaviour. 

These observations closely resemble the ideas developed to explain CdTe growth at a 



speed of ~ 0.5 ML/cycle in ||2^. Figure 4.12 shows Vg as a function of temperature. For 



350 K < T < 500 K, growth proceeds at a speed of about 0.5 ML/cycle. If temperature is 
increased further, Vg decreases quickly. Thus, our model is also capable of reproducing the 
experimental observation of plateaus in Vg{T). 

There are two differences between our computer simulations and experimental results. 
First, in experiments the plateau of ALE growth at 0.5 ML/cycle extends to considerably 
higher temperatures than in our model. Second, our model does not reproduce growth at a 
speed of one monolayer per cycle at low temperature. 

An inspection of the surface shows, that the breakdown of ALE at high temperature is 
mainly due to the quick evaporation of Te atoms during the Te phase. On the contrary, 
surfaces under a Cd flux appear to be comparatively stable. A simple consideration confirms 
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this picture. We expect a surface under a flux F of the particle species which terminates the 
reconstruction to be stable if the sublimation rate Tsub is smaller than F. If sublimation is an 
Arrhenius activated process with macroscopic activation energy -Esub ^-nd attempt frequency 
z^sub) tlie temperature Tjnst where rsub = F \s given by 

Tinst = - ,. f^]' y (4.11) 

A;ln(F/i/sub) 

Inserting the activation energy and the corresponding prefactor for the decay of a flat Te- 
terminated surface, we obtain Tjnst = 570 K which is close to the temperature 540 K at which 
Vg becomes zero. The experimental value of Tjnst = 700 K [^] corresponds to an activation 
energy of « 1.6 eVQ. Thus, the deviation from our value of 1.35 eV is about 20%. While 
this deviation might be cured by a tuning of model parameters, our considerations indicate 
an inconsistency between different experiments. Clearly, the fact that ALE is possible at 
temperatures above 600 K is incompatible with an activation energy as low as 0.96 eV for the 



decay of the Te terminated surface which has been measured in |11£] unless an unphysically 
low attempt frequency ~ 10^/s is assumed. 

The observation of a growth rate of IML/cycle has been explained with the formation 
of a surface reconstruction with a Te coverage of 3/2 in the Te phase. Here, one half of the 
potential Cd sites on the surface are occupied by Te* atoms. In the following Cd phase, 
these particles are bound to newly deposited Cd atoms where they provide sites which can be 
occupied by additional Cd atoms. Thus, in each cycle a complete CdTe monolayer is deposited 
in spite of the hard repulsion between Cd atoms on neighbouring sites in the y-direction. 

This effect can occur only if Te* atoms on a Te terminated surface are bound strong enough 
such that the rate of desorption of Te* is smaller than the particle flux. Since there is no 
interaction between Te* atoms in our model, macroscopical attempt frequency and activation 



energy are identical to the microscopical model parameters. Inserting these in equation 4.11 
we find that a layer of Te* atoms on the surface is stable at T < 116 K which is far below 
ALE growth temperatures. Simulations confirm this result. Conversely, the experimental 



ALE data and equation 4.11 yield an activation energy of 1.28 eV for the decay of a surface 
terminated by Te trimers. This high value indicates that there are attractive interactions 
between trimers. On the other hand, there must be also repulsive interactions which stop the 
Te coverage from becoming greater than 3/2. Thus, a consistent description of ALE growth at 
comparatively low temperature probably requires the introduction of anisotropic interactions 
between Te trimers. 

For completeness, we have also performed simulations of ALE on a stepped surface. At a 
terrace width of 16 lattice constants ALE growth proceeds in step flow mode. We find, that 
the growth rate Vg is identical to that found in ALE on a flat surface. However, there is no 
alternation between a rough and a flat surface in subsequent cycles. The deposited particles 
are captured by step edges and no islands are created. Therefore, we flnd a dominance of a 
c(2 X 2)cd arrangement in all Cd phases. 

4.9 Summary and discussion 

In summary, we have presented a model of growth and sublimation of CdTe(OOl) which 
might be an important step towards a realistic model of this surface. We have introduced an 

^ Since Tinst depends on ln(_F/i'sub), variations of particle flux and attempt frequency within a physically 
reasonable range have small effects. 
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efficient representation of the (001) surface of the zinc-blende lattice. On the one hand, the 
great advantage of a solid-on-solid model to allow for a two-dimensional representation of a 
three-dimensional crystal is preserved. On the other hand, it is possible to consider overhangs 
and atoms in the wrong sublattice as long as these are not incorporated into the crystal. The 
reconstructions of the Cd terminated surface are considered in a way which is inspired by 
the ideas developed in sections |2| and 0. Simulations of this model show, that a physically 
reasonable parameter set yields a rough quantitative agreement with experiments. 

To the best of our knowledge, we have performed the first simulations of atomic layer 
epitaxy. Our results show that it is possible to investigate the influence of surface reconstruc- 
tions on this important growth process by means of Monte Carlo simulations. Additionally, 
a comparison of the ALE growth rate as a function of temperature between simulations and 
experiments should be useful to fit model parameters. 

However, there are still some open questions. (1) What is the role of Te dimers on a 
non-equilibrium crystal surface? The lattice gas model of a flat CdTe(OOl) surface in thermal 
equilibrium (chapter ^) suggests that Te dimerization might be important, in particular at 
low pcd- As discussed in section 4/7, the dimerization of Te atoms might lower the Cd 
coverage of surfaces sublimating under vacuum at high temperature. This should lead to an 
increase of pcd and K'^^ and to a decrease of Kq^ under an external Cd flux. Experiments 
indicate that the Cd coverage under vacuum is indeed smaller than the values which we 
find in our simulations. In ||7^ pcd ~ 0.35 has been measured at the temperature of the 



crossover between dominant c(2 x 2)q^ and dominant (2 x l)cd- As discussed in section 2.3, 
this would correspond to a coexistence of a c(2 x 2)cd and a (2 x 1)tc phase on the surface 
in vacuum. (2) Our model does not reproduce ALE growth at a speed of one ML/cycle 
at low temperature. As discussed in section O, an understanding of this feature requires 
the consideration of interactions between Te trimers. Experiments have shown, that trimers 



arrange in rows parallel to the x-direction yielding a (2 x 1) reconstruction|21, p.l9| . There 
are no trimers on nearest neighbour sites in the y-direction which can be understood from 



the electron counting rule (section 1.1). This suggests a modelling by means of an attractive 
interaction between trimers on neighbouring sites in the x-direction and a hard repulsion 
between trimers on neighbouring sites in the y-direction. However, it is questionable whether 
the simulation of such a model is computationally feasible. Since Te* atoms are comparatively 
weakly bound, the diffusion of Te* atoms is fast. This leads to a dramatical slowing down 
of the simulations if there is a high density of Te* on the surface. Itoh et. al. [45, 46 1 faced 
a similar problem in the simulation of epitaxial growth of GaAs(OOl) where there is a high 
density of weakly bound As2 molecules. They solved it by treating all As2 states as one 
homogeneous particle reservoir without distinguishing between different arrangements of the 
molecules. Li our model, a similar treatment of the Te* atoms should be possible. The basic 
idea is to assume that the diffusion of Te* atoms is so fast that the probability of a Te* atom 
to interact with a given atom on the solid-on-solid surface does not depend on its history. 
Then, it is a reasonable approximation to assume that there is an equal density of Te* atoms 
above each site on the solid-on-solid surface. Interactions between Te* atoms might be treated 
in a mean field manner by assuming that the binding energy of a Te* atom depends on the 
total number of Te* atoms in the reservoir. 
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Chapter 5 



Singularity spectra and dynamic scale 
invariance of kinetically rough surfaces 



The past decade has raised much theoretical research in the subject of kinetic roughening of 
surfaces during growth. The investigation of this effect promises deep insight into statistical 
physics far from thermal equilibrium, see e.g. |12| for an overview. As discussed in section 



|1.5.1| , this process creates surfaces which have fractal properties on lengthscales smaller than 
the correlation length ^{t). 

We focus on a solid-on-solid model with a simple cubic lattice structure and nearest neigh- 
bour interactions only. The surface is represented by a two-dimensional array h of integers 
hx,y which denote the height of the column of atoms at {x, y). The rates of thermally activated 



processes are given by the Arrhenius law |1.2| . Energy barriers are determined by the number 
of bonds which must be broken to remove the particle from its initial site. 

Atoms may diffuse to nearest neighbour sites. The rate of such a process is i^exp(— (e^ + 
nen)/{kT)), where ej, and Sn are the binding energies of a particle to the substrate and to 
its n nearest neighbours in the same layer, u is the attempt frequency, and kT has its usual 
meaning. On each site, new particles are deposited with a rate F. There is no incorporation 
process. In contrast to earlier investigations of similar models [^], we permit the desorption 
of particles from the surface with rates z^exp(— (ei, -|- ne„)/(A;T)), where > Sb- Simulations 
of this model are performed by means of the continuous time algorithm which is introduced 
in appendix Clearly, such a simple model neglects many effects which are present in real 
materials. On the other hand, it captures the essential features of processes which occur on 
unreconstructed surfaces in the absence of a Schwoebel barrier. Therefore, it is adequate to 
obtain insight into fundamental properties of the roughening process which are not specific 
to a particular material. 

Recently, Arneodo et. al. [0, §, 2S, have suggested a method for the investigation of 
fractal and multifractal surfaces which uses the continuous wavelet transform. This is superior 
to the structure function approach which has to date solely been used in the analysis of models 
of epitaxial growth. In this chapter, we use the wavelet method for a precise characterization 
of kinetically rough surfaces. In particular, we investigate the influence of desorption on the 
surface morphology. We show that our model fulfils the dynamic scale invariance which was 



introduced in section 1.5.3 . We relate the dynamic exponents a, (3 with the static scaling 
properties of the multifractal surface. We conclude with some remarks on the question of 
universality of scaling exponents. 
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5.1 Anomalous scaling and multiscaling 



Originally, it was believed that the scaling properties of kinetically rough surfaces can be 
understood within the framework of Family- Vicsek scaling. Then, the properties of the sur- 
face are invariant under the dynamic scale transformation |1.9| . On lengthscales smaller than 
^{t) the surface is self-affine. Its Hurst exponent equals the dynamic exponent a. Under 
these conditions, the scaling properties of the surface are characterized completely by the two 
exponents a and (3 = a/z. These are assumed to be universal, i.e. independent of model 
parameters. 

Different methods have been suggested to measure a and j3 in computer simulations. 
Frequently, /5 is measured from the increase of the surface width with time, W{t) oc 



(equation |1.11 ). Since a characterizes both the static scaling properties of the surface at fixed 
time and the dynamic scaling properties, there are two fundamentally different approaches to 
determine its value. The static approach investigates the height-height correlation function 
G{1, t) at fixed time t. For small / := |/|, from equation 1.5 we expect it to increase ~ The 
dynamic^ approach analyzes the dependence of the surface width Wsat in the saturation regime 
on the system size N: According to equation 1.14, we have VFsat(-^) ~ An alternative 



method uses equation 1.12 : a and z are chosen such that the curves of G{1, t) versus l/t^^^ 



collapse on a unique function u within a large range of t and I. 

However, a careful analysis of simulation data 57 1 has shown that several models 

of epitaxial growth show significant deviations from this simple picture. First, frequently one 
obtains significantly different exponents from the static and the dynamic approach. This 
phenomenon is commonly denoted as anomalous scaling. Second, a powerlaw behaviour of 
G{l,t) for small |/| does not necessarily imply that the surface is self-affine. The standard 
approach to test for self-affinity considers the height-height correlation functions of order q 



T{q,l,t) 



f{x,t)-f{x + l,t) 



(5.1) 



which are generalizations of G{l,t) = T{2,l,t). Here, f{x,t) is the reduced surface height 
h{x,t) — {h{x,t)). On a self-affine surface we expect T{q,l,t) ~ P'''^ for small I, where 7^ = a 
for all q. However, sometimes one encounters the phenomenon of multiscaling: the initial 
powerlaw behaviour of T{q, I, t) yields a hierarchy of g-dependent exponents 7^. 

These observations can be interpreted within the mathematical framework of multifrac- 
tality: The Holder exponent |l^, |l^, 72 1 xi^o) of a function / at xq is defined as the largest 
exponent such that there exist a polynomial of order n < xi^o) and a constant C which 



yield \ f{x) - P„(f - xo)| < C\x - xq 



in the neighbourhood of xq. The Holder exponent 



characterizes local symmetry properties of the surface. In the vicinity of xq, we have 

f{xo + bl, t) - f{xo, t) ~ b^^^"^ [f{xo + - f{xo, t) 

Thus, the Holder exponent is a local counterpart of the Hurst exponent: a self-affine function 
with Hurst exponent H has xi^) = H everywhere. However, in the case of a multiaffine 
function different points x might be characterized by different Holder exponents. In this case, 
the surface consists of interwoven fractal sets of points where there is one particular Holder 
exponent. This general case is characterized by the singularity spectrum -D(x)i which denotes 
the Hausdorff dimension of the set of points where x is the Holder exponent of /. 

^In the literature (e.g. | ]6o| , ^) the static and the dynamic approach are frequently denoted as 

"local" and "global" approach, respectively. 
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5.2 The wavelet approach to multifractality 



There is a deep analogy between multifractality and thermodynamics |T2| , 105 , 125 |, where 
the Holder exponents play the role of energy, the singularity spectrum corresponds to entropy, 
and q plays the role of inverse temperature. So, theoretically D{x) might be calculated via a 
Legendre transform of 7^: D{x) = mmq{qx — qjq + 2) |7^, a method which has been 

called structure function (SF) approach. However, its practical application raises fundamental 
difficulties: First, to obtain the complete singularity spectrum, one needs 7^ for positive and 
negative q. But as \ f{x,t) — f{x + l,t)\ might become zero, T{q,x,t) is in principle undefined 
for q < 0. Therefore, only the left, ascending part of D{x) is accessible to this method. 
Additionally, the results of the SF method can easily be corrupted by polynomial trends 
in f{x) It might be due to these difficulties, that — to our knowledge — no attempt to 

determine the singularity spectrum of growing surfaces from 7^ has ever been made. On the 
basis of simulations of several simple growth models, it has been argued that the 7^ might 
collapse onto a single 7 in the limit t — > 00, which characterizes the asymptotic universality 
class of the model This would imply that multiscaling is a finite size effect. 

However, there is no physical reason why this should be the case. A precise measurement of 
D(x) will help to test this hypothesis. 

To this end, we follow the strategy suggested by Arneodo et. al. [Q, IC, 121 1, which circum- 
vents the problems of the SF approach and permits a reliable measurement of the complete 
singularity spectrum. Mathematically, the wavelet transform of a function f{x) of two vari- 
ables is defined as its convolution with the complex conjugate of the wavelet il^i^x) which is 



dilated with the scale a and rotated by an angle 9 |121]: 



T4f]{b,e,a) = C-'^'a--' I d'x r{a^''i^-e{x - b))f{x). 



(5.2) 



Here is the usual 2-dimensional rotation matrix, and = (27r)^ J d'^k\k\~'^\'ijj{k)\'^ is 
a normalization constant, whose existence requires square integrability of the wavelet ^lJ{x) 
in fourier space. Apart from this constraint, the wavelet can (in principle) be an arbitrary 
complex- valued function. Introducing the wavelet ips{x) = S{x) — S{x + n), where n is an 
arbitrary unit vector, one obtains easily 



C, 



-1/2 



dH 



f{b)- f{b + aR.gn) 
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Tji,s[f]{b,6,a) (xr{q, alien) 



(5.3) 



Consequently, a calculation of the moments of the wavelet transform of the surface yields 
the SF approach as a special case. To avoid its weaknesses, two major improvements are 
necessary. 

First, we use a class of wavelets with a greater number of vanishing moments than 
ips{x)- Then, the wavelet is orthogonal to any polynomial of degree less than n^. This 
increases the range of accessible Holder exponents and improves the insensitivity to polynomial 
trends in f{x). We introduce a two-component version of the wavelet transform 




(5.4) 
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where the components ^i, ^'2 of the analyzing wavelet ^ are defined as partial derivatives 
of a radially symmetrical convolution function $(x): ^i{x) = d^/dx, ^'2(2?) = d^/dy. Then 
T^[/](6, a) can be written as the gradient of /(x), smoothed with a filter $ with respect 
to b. This definition becomes a special case of equation |5.2| , when multiplied with fig = 
{cos{9),sm{6))~^ , yet allows for an easier numerical computation^. For example, <1> can be a 
gaussian, where = 1, or $i(x) = (2 — x^) exp(— which has two vanishing moments. 

Second, the integration over b in equation |5.3| is undefined for q < 0, since the wavelet 
coefficients might become zero. The basic idea is to replace it with a discrete summation over 
an appropriate partition of the wavelet transform which takes on nonzero values only, but 
preserves the relevant information on the Holder regularity of /(x). In the following, we give 
a brief outline of the rather involved algorithm and refer the reader to [0, IC , 121 ] for more 
details and a mathematical proof. 

The wavelet transform modulus maxima (WTMM) are defined as local maxima of the 
modulus M^[f]{b,a) := \ f^{f](b,a)\ in the direction of f^[f](b,a) for fixed a. If <I>(x) is the 
gaussian function, the computation of their positions is equivalent to the detection of edges in 
image processing ||6^, Similarly, for arbitrary wavelets the WTMM lie on connected curves 
which trace structures of size ~ a on the surface. The strength of each curve is characterized 
by the maximal value of M^[f](b,a) along the line, the so-called wavelet transform modulus 
maxima maximum (WTMMM) 0. While proceeding from large to small a, successively 
smaller structures are resolved. Connecting the WTMMM at different scales yields the set C 
of maxima lines I, which lead to the locations of the singularities of f{x) in the limit a — > 0. 
The partition functions 

Z{q,a)= ( ^ s^iP M^[f]{b,a')\ ~ a^^^) for a^O (5.5) 

leC{a) \(b,a')el,a'<a ) 

are defined on the subset C-ia) of lines which cross the scale a. From the analogy between 
the multifractal formalism and thermodynamics, -D(x) is calculated via a Legendre transform 
of the exponents r(g) which characterize the scaling behaviour of Zi^q^ a) on small scales a: 
D{x) = mi%('?X ~ '^{q))- Additionally, T{q) itself has a physical meaning for some q: — r(0) 
is the fractal dimension of the set of points where h{x) < 00, while the fractal dimension of 
the surface f{x) itself equals max(2, 1 — t(1)). 

5.3 Results 

In our simulations, we choose the parameters u = 10^^ s, £{, = 0.9 eV, e„ = 0.25 eV, and a 
temperature T = 450 K. To study the influence of desorption, we consider three models with 
different activation energies e^: in model A desorption is forbidden, i.e. = 00. Models B 
and C have = 1. 1 eV and Sy = l.O eV. We simulate the deposition of 2 • 10^ monolayers 
at a growth rate of one monolayer per second on a lattice N x N unit cells using periodic 
boundary conditions, our standard value being = 512. To check for finite size effects, we 
have also simulated N = 256. In all presented results averages over 6 independent simulation 
runs have been performed. Snaphots of the surface in a simulation of model A are shown in 



figure 5.1. Visually, surfaces of models B and C look quite similar. 



^For simplicity, the irrelevant constant has been omitted. 
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t = 200s 



t = 2000s 



t = 20000s 



Figure 5.1: Evolution of the surface in a simulation run of model A. Brighter shades of grey 
denote greater surface heights. 
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Figure 5.2: Panel (a): r(g) in model A as obtained from investigations with wavelets derived 
from the filters (+), <l>i (x), $2 (Q) and $3 (■). Note the deviations of the data obtained 
with the gaussian function ^q. Panel (b): r(g) for models A (+), B(x) and C (□). All data 



have been obtained from surfaces after 2 
the order of symbol sizes. 



10^ s of simulated time. Sizes of errorbars are on 



5.3.1 Singularity spectra 

First, we have checked our results for artifacts resulting from properties of the analyzing 
wavelet rather than the analyzed surface by using different convolution functions <^„: 
is the gaussian function, ^n,n > 1 are products of gaussians and polynomials which have 
been chosen in a way that the first n moments vanish. Then, the analyzing wavelets have 
= n + 1 vanishing moments. We find (figure ^.2^ ), that the T{q)-cmve obtained with 
$0 deviates significantly from those obtained with $1, $2 and $3. The latter agree apart 
from small differences which are mainly due to the discrete sampling of the wavelet in the 
numerical implementation of the algorithm. This is explained by the theoretical result ||7^ 
that dT{q)/dq = for q < ^crit. < if the number of vanishing moments of the analyzing 
wavelet is too small. Consequently, the agreement of the other curves proves their physical 
relevance. 



Figure 5.2b shows averages of T{q) curves obtained with the convolution functions ^i, <I>2 
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$3 from surfaces after 2 • 10^ s of growth on an initially flat substrate. For all our models, 
their nonlinear behaviour reflects the multiaffine surface morphology. From the fact that these 
curves are reproduced within statistical errors in simulations with = 256 we conclude that 
finite size effects can be neglected. Clearly, desorption reduces the slope of r(g), although only 
a small fraction of the incoming particles is desorbed: 0.18% in model B and 2.57% in model C 
with slightly higher values at earlier times. The corresponding singularity spectra are shown 
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Figure 5.3: Panel (a): Singularity spectra obtained from a Legendre transform of the data 



in figure 5.2b. Panel (b): Data collapse of rescaled PDFs of surface heights for times between 
150 s and 2 • 10^ s (model A) respectively 150 s and 7500 s (model C) . The corresponding scaling 
exponents are summarized in table |5.1| . Time is measured in seconds. In both panels, the 
symbol + represents results of model A, x corresponds to model B and □ denotes results of 
model C. Model B which is not shown in panel (b) yields a data collapse of similar quality as 
model A. 



in figure 5.3a. They have a typical shape whose descending part seems to be symmetrical to 
the ascending part and which changes at most slightly, while the whole spectra are shifted 
towards smaller Holder exponents as desorption becomes more important. We emphasize 
that we find no evidence for a time dependence of the singularity spectra within the range 
9700 s < t < 2 - 10^ s, such that our results do not support the idea of an asymptotic regime 
which is governed by Family- Vicsek scaling. However, the accessible time range of computer 
simulations is limited, so we cannot finally disprove the existence of such a regime. 



5.3.2 The relation between static and dynamic scaling properties 

In Family- Vicsek scaling, we have H = a such that the scaling exponent of a static surface 
is identical to one of the dynamic exponents. The multifractal formalism has replaced the 
unique Hurst exponent H of a self-affine surface with a wide spectrum of Holder exponents. 
This raises the question whether the analogy between static and dynamic scaling behaviour 
holds also in the presence of multiscaling such that there is a distribution of dynamic scaling 
exponents. Alternatively, it is possible that multiscaling affects only the properties of static 
surfaces while dynamic scaling is governed by two unique exponents a and /3. 

To answer this question, we need a reliable method to decide whether the properties of the 



surface are invariant under the transformation L£. To this end, we investigate the probability 
distribution function (PDF) Pf{f,t) of the reduced surface height. Invariance under dynamic 
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Figure 5.4: Data collapse of -Pa/- Panel (a) shows data from model A in the time range 
150s < t < 2 • 10^ s at Ijt^!'' = 0.33 (+) and //t^/^ = 8.42 (©). In panel (b), we show results 
for model C within the time range 150s < t < 7500s at //t^/^ = 0.63 (+) and l/t^l^ = 15.1 
(O). The data collapse is good apart from small lattice effects at A/ = 0. Model B yields 
results similar to model A. In all results shown, the vector / is parallel to the x-direction. 



scaling implies that Pf{b°'f,bH)d{b°'f) = Pf{f,t)df. This is fulfilled if 

^/(/>0 = ^P/(^)- (5.6) 

Consequently, within a large time range the rescaled PDFs Pjt^ should collapse onto a single 
function P/ if they are plotted as a function of f /t^. Additionally, we consider the PDF 
P^f{Af,l,t) of the height increment at a distance /, A/ := f{x,t) — f{x + l,t). Similar to 
equation we obtain 

^A/(A/,M) = ^PA/(^^,^^. (5.7) 

Therefore, in a large time range one should obtain a data collapse if t^PAf is plotted versus 
Af /t^ and the ratio l/t^/^ is kept fixed. Surface width and height-height correlation functions 
of arbitrary order are moments of Pj and Pa/? respectively. Consequently, our analysis of the 
PDFs considers all the information which is contained in W{t) and T{q,l,t) for arbitrary q. 

To obtain values of the dynamic exponents a and (3, we apply standard methods which 
have been reported in the literature. We measure (3 from the increase of the surface width with 
time, which follows a power law for t > 150 s in models A and B and for 150 s < t < 7500 s in 
model C which then starts to approach the final saturation regime. The dynamic exponents 
a and z have been obtained from the data collapse of the scaled height-height correlation 



function G{l,t). The results are summarized in table 5T 



In figure |5.3[ b, we show the data collapse of Pj in models A and C which is of high quality. 



Model B yields similar results. Examples of rescaled Pa/ are shown in figure The data 
collapse is good for A/ ^ 0. However, the probability that the height difference between 
points at a distance l/t^^^ is zero does not follow the behaviour which is expected from the 



76 



assumption of dynamic scale invariance. This deviation is a finite size effect. Dynamic scale 
invariance is a continuous symmetry property. On the other hand, due to the discrete lattice 
structure, differences between surface heights can assume only integer values. This makes it 



impossible to fulfil dynamic scale invariance exactly, since the transformation |l.9| maps integer 
height differences to noninteger height differences. If we approximate a real random number 
r with an integer i, on average the relative error \i — r\/r is greater for smaller r. Therefore, 
it is plausible that the deviations of the empirical P/\f{Af,l,t) from the theoretical result in 
equation are greater for small A/. However, in the limit of late time and lengthscales 
which are large compared to the lattice constant, corrections due to the lattice structure can 
be neglected. Consequently, there is no indication that there is a violation of dynamic scale 
invariance apart from small finite size corrections. This parallels the result of Krug [^] that 
there is an unique exponent /? in the one-dimensional Das Sarma-Tamborenea model. 

Finally, we relate the static scaling properties of the surface which have been determined 
precisely by means of the WTMM method with the dynamic scaling properties. Our results 
(table |5.l| ) show, that the dynamic scaling exponent a and the value Xm of the Holder exponent 
which maximizes D{x) agree within errorbars. This empirical finding can be explained with 
a saddle point argument. We calculate the standard deviation of the surface in a square of 
size L X L, 



wiL, tf = £ <fxf{x, t? = ^ j dx £ ^ 



L 

(fx5{x-x{S))f{x,tf . (5.8) 



Kx) 

Since I{x) grows like L^^^'^ with L, for large L the integral over x will be dominated by I{xm)- 
Thus, w{L, t) is governed by the subset of points which has the greatest fractal dimension. 
Consequently, in a multiaffine surface w{L,t) is identical to that of a self-affine surface with 
Hurst exponent Xm- In particular, for L < ^{t) we have w{L,t) ^ L^™ (equation |1.7D . On 
the other hand, dynamic scale invariance implies w{L,t) ~ for L < ^(t). Consequently, 
Xm = a- 

The conventional picture of anomalous scaling |6C, 62, B^] notes the difference between a 



"global a" (the dynamic exponent) and a "local a" which is determined from the power-law 
behaviour of G{l,t) for small I and, in our notation, is termed 72. Within the multifractal 
formalism, the local a simply corresponds to some Holder exponent on the ascending part of 
the singularity spectrum. In contrast to the dynamic scaling exponents, it is not related to 
any symmetry property of the model itself. Therefore, we expect the dynamic exponents a 
and z to be most significant to characterize the scaling properties of a growth model. 



5.4 Conclusions 

Our results support the following picture of the scaling behaviour of kinetically roughening 
surfaces: Dynamic scale invariance is fulfilled on time- and lengthscales where lattice effects 
can be neglected. Surfaces at fixed time have multifractal properties which are characterized 
by a wide spectrum of Holder exponents. The Holder exponent which maximizes D{x) is 
identical to the dynamic exponent a. There is no indication that the singularity spectra 
narrow with time. Consequently, it is plausible that the multifractal properties of the surface 
are preserved in the limit of infinite time and there is no transition to Family- Vicsek scaling. 
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Model 


Pd 


/3 


Xm 


a 


z 




A 





0.19 ± 0.01 


0.54 ± 0.01 


0.55 


2.9 


2.32 ± 0.01 


B 


0.18 % 


0.17 ± 0.01 


0.52 ± 0.01 


0.51 


3.3 


2.35 ± 0.01 


C 


2.57 % 


0.11 ± 0.01 


0.38 ± 0.01 


0.39 


3.5 


2.45 ± 0.02 



Table 5.1: Simulation results: pd is the fraction of particles which desorbs, Xm the Holder 
exponent which maximizes D(x), (3 and z are the dynamic scaling exponents, and Df 
the fractal dimension of the surface. 



Our discussion of equation 5.5 indicates that this is the most general scaling behaviour 
given the following conditions are fulfilled: (1) The statistical properties of the model are 
invariant under the dynamic scale transformation (2) The surface is scale-free such that 
the statistical properties of sections of the surface of size L are independent of L whenever 
L < £,{t). Then, ^(t) is the only relevant lengthscale of the surface. This picture contains 
Family- Vicsek scaling as the special case where there is only one Holder exponent. 

According to Lopez |^], there are continuum models of epitaxial growth which exhibit 
anomalous scaling but not multiscaling. In particular, this is the case for the Mullins equation 
and the Lai-Das Sarma equation. Since the singularity spectrum is the Legendre transform 
of 7g, in this case D{x) has a single peak at one characteristic Holder exponent which equals 
the exponent 72. Thus, the static scaling properties which are measured on short lengthscales 
(equation 5.5) are not related to the dynamic scaling properties. This implies necessarily, 
that the dynamic properties are determined by large structures which appear smooth on 
small lengthscales. Consequently, the surface is not scale-free. Instead, we expect that there 
is a second lengthscale E{t) < ^{t). On lengthscales ^ S(i)> the surface appears self-affine 
with a Hurst exponent = 72. The dynamic properties are determined by structures on 
lengthscales between and ^(t). Note, that dynamic scale invariance of the model implies 

Table |1 
results in IBS 



summarizes our quantitative results. Model A without desorption reviews the 
, which have been obtained with slightly different activation energies on smaller 
systems and shorter timescales. Models B and C show, that desorption is an important 
process, which, although it affects only a small fraction of the adsorbed particles, must not be 
neglected, since it alters the scaling properties of the surfaces by reducing both a and (3. Since 
the dynamic exponents depend continuously on the desorption rate, the paradigm of a few 
universality classes characterized by a small number of exponents which are independent of 
details of the model is not adequate to catch the features of kinetic roughening in our computer 
simulations. Since the time range of computer simulations is limited by the available CPU 
time, our results cannot finally disprove that there is a crossover to a different, possibly 
universal scaling behaviour after the deposition of a number of monolayers which is large 
compared to our 2 • 10^. However, if this should indeed be the case, it is improbable that the 
asymptotic scaling behaviour might ever be observed on experimentally relevant timescales of 
a few hours of growth. 
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Chapter 6 

The influence of the crystal lattice on 
coarsening in unstable epitaxial growth 



In spite of considerable efforts |48, 69, 95, 110, 111, 117, 118 |, see e.g. [S6| for an overview, 
a thorough theoretical understanding of the coarsening process of mounded surfaces is still 
lacking. In this chapter we investigate the problem of growth in the presence of a strong 
Schwoebel barrier which hinders interlayer transport. As discussed in section 1.5.2 , this leads 
to an instability of the flat surface. During an initial transient, mounds form on the surface 
which start to merge as soon as their flanks have assumed the stable slope. It is generally 
accepted that this asymptotic coarsening regime fulfils the dynamic scale symmetry which 
has been introduced in section 1.5.3. 



It was first pointed out by Siegert et. al. |110| , 111| that lattice symmetries may play an 
important role in the coarsening process. These authors investigated continuum equations, 
using an analogy between coarsening and a phase ordering process which has recently gained 
popularity ||5^: areas of constant slope should correspond to domains of a constant order 
parameter. They derived scaling exponents a = l,/3 = 1/3 on surfaces with a triangular 
symmetry, and (3 = l/(3\/2) ~ 0.24 for generic cubic surfaces, while [3 = 1/3 requires a 
fine-tuning of parameters. Moldovan and Golubovic applied similar methods to surfaces 
with hexagonal symmetry and found (3 = 1/3. 

However, Monte Carlo simulations have raised doubts on these predictions, since they 
yield /3 ~ 1/3 on cubic surfaces for a great range of parameters, see e.g. ||6|, 101, 104 1. In 
this chapter, we adress the following questions: (1) What are the mesoscopic processes which 
make the mounds coarse? (2) How does the coarsening process depend on the crystal lattice 
and its symmetries? (3) Do these results support a deeper analogy between coarsening and 
phase ordering? 



6.1 Effective single particle simulations 

To answer these questions, we perform computer simulations of growth on (001) surfaces of the 
simple cubic (sc) lattice, the simple hexagonal lattice (sh), the body-centred cubic (bcc) lattice 
and the (0001) surface of the hexagonal close packing (hep). This is done under solid-on-solid 
conditions, i.e. the effects of overhangs or dislocations are neglected. Then, the simple lattices 
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can be represented by a square (sc) and a triangular mesh (sh) of integers respectively, which 
denote the height h{x) of the surface. We build the bcc (hep) lattice out of two intersecting sc 
(sh) sublattices, which contain the even and the odd heights, respectively]^. Here, an adatom in 
a stable configuration is bound to 4 (3) neighbours below in the other sublattice, which will be 
denoted as vertical neighbours in the following. Particles with fewer vertical neighbours form 
overhangs which are forbidden by the solid-on-solid condition. This is physically reasonable, 
since such particles are only weakly bound and therefore these configurations will be unstable. 
Additionally, particles may have neighbours in the lateral direction which are in the same 
sublattice. 

The investigation of the coarsening process requires a fast algorithm which allows for the 
simulation of the deposition of thick films on comparatively large systems. Kinetical Monte 
Carlo techniques which consider the moves of many particles on the surface simultaneously 
are too slow for our purpose even if the efficient continuous time algorithms introduced in 
appendix M are used. Instead, we simulate the moves of a single particle from deposition until 



an immobile state is reached [15, BQ, 101]. 



An adatom impinges on a randomly chosen lattice site. As discussed in section 1.4.1, it is 



reasonable to assume that a newly deposited adatom undergoes an incorporation process. In 



our model, the particle funnels downhill p2l , |33| , 13C ] to the lowest (vertical) neighbour site. 
On bcc and hep this is repeated until it reaches a stable site as defined above. On the simple 
lattices this is a site x where all nearest neighbour sites have a height > h{x). 

Then, the adatom diffuses on the surface. If the particle has no lateral neighbours, one of 
its neighbour sites is chosen at random. On the simple lattices, the particle moves to this site 
only if its height does not change, i.e. we introduce an infinite Schwoebel barrier. On the bcc 
or hep lattice, the particle is moved to the neighbour site if it is stable. This condition implies 
an infinite Schwoebel barrier, too. Similar to our model of the zinc-blende lattice introduced 
in chapter ^ interlayer diffusion would require the consideration of additional hopping vectors. 

As discussed in section L5, from time to time adatoms which diffuse around on a flat 
surface collide and nucleate a new island. The typical distance li between neighbouring 
islands depends on the ratio of the adatom diffusion constant D and the particle flux F. In 
our effective single particle dynamics, we consider this effect by limiting the maximal number 
of diffusion steps of a particle. If an atom has performed l'^ diffusion steps, diffusion stops and 
the particle becomes an island nucleus which will capture subsequent particles. The dijfusion 
length Id fixes the number of islands in the first layer which is identical to the initial number 
of mounds. In the later stages of growth, the typical terrace width is much smaller than 1^ 
and nucleation occurs only at the top terraces of the mounds. 

As soon as a particle has a lateral neighbour, it is bound to it. This process is irreversible 
in the sense that we forbid diffusion processes which reduce the number of bonds like the 
detachment from a step edge. In section L5 we have shown that this is a reasonable ap- 
proximation if the rates of these processes are much smaller than F. However, the adatom 
may diffuse along the edge. After steps, or if the particle has reached a kink site, it is 
fixed to the surface. If not stated otherwise, diffusion around corners is allowed. Analogous 
to planar diffusion, Ik determines the typical distance of nucleation events in the effectively 
one-dimensional step edge diffusion, see [124] for a discussion. 

In this model, we measure time t in units of the time needed to deposit one monolayer 

^For simplicity, we assume the spacing between the layers to be one unit length. Our algorithm depends 
only on the topology of the lattice. 
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200 ML 



2000 ML 



t = 20000 ML 



Figure 6.1: Coarsening of mounds on the hcp(OOOl) surface at M = 300, Id 
Brighter shades of grey correspond to greater values of the surface height. 



15 and L = 20. 



(ML). This algorithm may be programmed very efficiently and is at least one order of magni- 
tude faster than full diffusion Monte Carlo algorithms. 

The two essential simplifications in our model are (a) the effective representation of nucle- 
ation events in a single particle dynamics and (b) the consideration of irreversible processes, 
i.e. infinite energy barriers for detachment and downward diffusion at edges. For the simple 
cubic lattice, the model was studied in [15, 101, 104| with particular emphasis on the role of 



step edge diffusion. In | 104 | a full diffusion model is considered for comparison which includes 
detachment as well as finite Schwoebel barriers and displays the same asymptotic coarsening 
exponents as measured in the simplified model. 

In all our simulations we choose Id = 15 which controls the initial island distance. We 
simulate a variety of different values of the step edge diffusion length in the range between 
Ik = 1 and Ik = 20. The simulations are performed on a square oi N x N lattice constants 
using periodic boundary conditions (bcc and sc) and a regular hexagon with edges of length 
M (hep and sh) using helical boundary conditions, our standard values being = 512 and 
M = 300. For every parameter set we have performed 7 independent simulation runs. 



6.2 Surface morphologies 

During the first (about 100) monolayers of growth on an initially flat surface islands nucleate 
on which mounds build up and take on their stable slope. Then, the asymptotical coarsening 



regime starts. As an example, in figure 6.1 we show the evolution of the hep surface at Ik = 20. 
As expected, on the simple lattices the mounds obtain regular shapes which are determined 
by the symmetry of the surface: square pyramids on sc, hexagonal ones on sh (figure |6.2|) . 
On bcc and hep however, rounded corners as well as sharp corners can be found (figures |6.1| , 
and S.5). 



Here one finds two fundamentally different types of mounds: On the one hand there are 
mounds, where all corners are rounded with octagonal shapes on bcc (figure |6.5| b), and 12- 
cornered ones on hep. On the other hand, one observes a breaking of the symmetry given by the 



lattice: approximately triangular mounds on hep (figure 3.4 b) and oval shapes on bcc (figure 
|6.5| a), where every second corner is sharp while the others are rounded. We have observed 
that in the process of coarsening, on all lattice structures mounds merge preferentially at 
the corners, while mounds touching at the flanks are a metastable configuration. On bcc 
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[100] 



Figure 6.2: Contour plots of surfaces of the simple lattices at t = 2 • 10^ ML. Brighter 
contours denote greater surface heights. The arrows show the average diffusion current on 
the surface during the deposition of additional 200 ML. The figures show sections of 250 x 250 
lattice constants of systems of size N = 512 (simple cubic, panel (a)) and M = 300 (simple 
hexagonal, panel (b)). Diffusion length and step edge diffusion length are Id = 15 and If^ = 20, 
respectively. 
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Figure 6.3: Histogram of the slopes on 
simple cubic surfaces (iV = 512, la = 15, 
If^ = 20, t = 2 • lO'' ML). High probabili- 
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Figure 6.4: Simulation results on hep (M = 300, Id = 15, h = 10, t = 20000 ML). Panel (a): 
Histogram of slopes (density plot). High probabilities are drawn dark. Panel (b): Contour 
plot of a part of size 250 x 250 lattice constants of the surface. High levels are plotted in light 
grey. The arrows show the diffusion current as measured during the deposition of additional 
200 ML. 



and hep, the metastable configuration consists of mounds with completely rounded corners. 
Mounds whose shapes break the lattice symmetry are always in the course of merging with a 
neighbour they touch with a sharp corner. Note that similar configurations were observed in 
a full diffusion simulation of bcc(OOl) growth [^]. 

To characterize the morphology quantitatively, we investigate the slopes m = V/i of the 
surface. On average, |m| obtains a value which can be estimated by a simple argument |84]: 
The Schwoebel effect yields an uphill current due to the preferential attachment of particles 
to steps from the lower terrace which is compensated by a downhill current due to tunneling. 
If we assume that a particle reaches its final height after at most one incorporation step, 
the sc and the sh surface select \fh\ = 1/2, while bcc and hep select |m| = 1. In practice, 
the selected slopes are slightly smaller, since the maximal number of incorporation steps in 
the funneling process is unlimited. Since a direct computation of the numerical gradient 
of simulated surfaces fails due to the discrete heights, we first apply a gaussian filter with 
variance a = A. This value has turned out to be a good compromise: it removes the atomistic 
structure but preserves the shape of the mounds. On sc and sh, we find a regular square and 
hexagon of slightly broadened peaks, respectively, which correspond to the slope of the flanks 
of the mounds (figure |6^ ). This is what one expects at pyramidal mounds. 

On bcc and hep, however, we obtain a more complicated pattern. Figure 6^a shows a 
histogram of the slopes on hep surfaces: due to the round shapes of the mounds, one finds a 
pronounced maximum of the probability density function in every direction of rh. The most 
probable value of |m| has a directional dependence: it is minimal in the lattice directions 
(0.88 ± 0.02), and maximal (0.99 it 0.02) in the intermediary directions. Consequently, the 
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Figure 6.5: Contour plots of parts of size 250 x 250 lattice constants of bcc surfaces (A^ = 512, 
Id = 15, Ik = 10, t = 2 • 10^ ML). Panel (a): section dominated by quickly coarsening mounds 
breaking the surface symmetry. Panel (b): metastable configuration. 



rounded corners are steeper than the flanks. This is possible because the rounding occurs at 
the base of the mounds, whereas close to the tips corners are typically sharp, cf. figure 6.5 . 
On bcc the results are equivalent, apart from the different symmetry. 



6.3 Difrusion currents 

In order to gain deeper insight into the coarsening process, we investigate the material trans- 
port on mesoscopic lengthscales by tracing the motion of particles on the surface. On all 
lattice structures, we observe a strong uphill current at the corners of mounds which is a 
geometrical effect: an adatom attaching at the corner is transported over a mean distance l^, 
namely with equal probability along either of the two flanks of the mound. This results in 
a net current j in the direction of the bisector of the angle at the corner, which is directed 
uphill. A simple calculation yields 

IjI = rah cos - 

if there is a straight step edge of length ^ Ik, which is the case on the simple lattices (figure 
[6.21) . Here, is the rate at which particles attach to a step edge, cf) is the angle between the 
two steps which meet at a corner. We have (p = Tr/2 on surfaces with cubic symmetry and 
(j) = 27r/3 on hexagonal surfaces. 

On the bcc and the hep lattice however, material transport along step edges is severely 
restricted due to the rounded corners. There, particles diffusing along the edge are typically 
reflected. This can be understood from the fact that the slope of the rounded corners is high. 
Consequently, there are sites next to a step where an additional atom would have less than 
4 (on bcc) or 3 (on hep) vertical neighbours. Since we consider only diffusion to stable sites, 
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an atom cannot cross such sites. 

In the metastable configuration, this leads to a current towards the middle of the smooth 
step edges at the flanks of the mounds which have an extension ~ (figure S.5). This current 
has a non- vanishing uphill component due to the Schwoebel effect. On the rounded corners 
on the contrary, one observes a weak downhill current. This is explained by the observation 
of steep gradients (figure |6.4| a) , at which the downhill current due to funneling dominates the 
uphill current of the Schwoebel effect. On these surfaces the stable slope at which the currents 
cancel is not assumed locally, but only in the global average, which results in spatial current 
patterns on the surface. An analogous behaviour is observed at the symmetry-breaking. 



merging mounds (figures 6.4b, |6.5| a). Here, the smooth edges are in the proximity of the 
sharp corners. Particles are reflected at the rounded corners, too. This yields a net current 
towards the sharp corners which compensates the geometrically induced uphill current exactly 
if the smooth edges have a length Z^. This fact explains the broken symmetry of these mounds: 
if there is any process which transports matter towards sharp corners, this material is not 
completely diffusing inward as it would be the case on pyramidal mounds but is attached 
to the corner and makes it overgrow its neighbours. Merging of mounds, however, is such a 
process: consider a terrace which goes around two mounds touching each other at the corners. 
Due to its curvature, in this contact zone there is a high concentration of kink sites which 
make it a sink for diffusing particles. 



6.4 Theoretical estimation of scaling exponents 

Our observation that mounds merge perferentially in positions with touching corners gives 
strong evidence that the current which fills the gap between two mounds plays a dominant 
role in the process of coarsening, at least in the case of large Ik where it is efficient. If this is 
the case, then a simple consideration yields the scaling exponents: Due to step edge diffusion, 
adatoms are transported over a typical distance Ik- If h is noticeably greater than a few 
lattice constants, this is much more than the average terrace width. Then, material transport 
via diffusion on terraces is small compared to transport via step edge diffusion and can be 
neglected. Our investigation of the diffusion currents shows, that the current into the gap is 
significant on a few atomic layers in the vicinity of the contact point only. In consequence, this 
current is independent of the size of the mounds if the latter is much greater than Ik- Since 
the volume of the gap is proportional to L^, if L is the typical distance between the mounds, 
it will take a time t ~ to fill it. In consequence, L t^^^ =^ z = 3. Since a = 1 in the 



presence of slope selection, this yields /3 = 1/3. In contrast to the effects discussed in [110, 111] 



these considerations are completely independent of lattice symmetries. In [117, 118[ a similar 



argument has been applied to the case of coarsening in the absence of dominant step edge 
diffusion processes ("bond energy driven coarsening"), where material is transported mainly 
by terrace diffusion. This yields (3 = 1/4. The same exponent is obtained, if coarsening is 
exclusively due to fluctuations in the particle beam ("noise assisted coarsening"). We expect, 
that these processes dominate step edge diffusion in the limit of small Ik which leads to a 
transient from f3 = 1/4 to f3 = 1/3 if /yt is increased. 
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6.5 Measurement of scaling exponents 



We apply a variety of methods to measure the scahng exponents. This is important as 
a consistency check and to ehminate systematical errors which might be intrinsic to some 



methods. According to equation |1.11| , P can be obtained from the increase of the surface 
width W{t) ~ t^^ PI. This power law behaviour may be corrupted by the presence of noise 
on the surface profile, an additional contribution to W{t) which has been called intrinsic 
width 1 12, E^. Similarly, the number of mounds in the system will decrease like if 



this scaling hypothesis holds. We measure n.^ by counting the number of top terraces in the 
system. Since a single particle is counted as a top terrace, this method may be misleading if 
particles nucleate on terraces at the flanks of mounds. A method which is widely used in the 
literature uses the fourier transform f{k) of the reduced surface height f{x): The structure 
factor S{k) := (^f (k) f (—k)"^ is maximal at nonzero wavenumbers \km\ = 2-n /Im, if there are 

structures at a typical lengthscale Im on the surface. Since Im ~ t^^^ , \km\ will decay ~ t~^/^ . 
In practice, a direct search for the maximum often fails due to noise effects; one avoids this 
problem by calculating the averages kt^ := {Y.tMS{kY^ /J2kS{ky. However, the choice 



of the correct power p is a bit arbitrary. 
6.5.1 The wavelet method 

To avoid these difficulties, we propose a method which uses the continuous wavelet transform 
which was introduced in chapter ||. Due to the easier numerical computation, we employ the 
two-component version ^.4| . 

The basic idea of this transformation is to convolve h{x) with the wavelet which is dilated 
with the scale a. As already noted in section ^.2[ in a wavelet transform where the radially 
symmetrical filter <l>(x) is a gaussian, the search for the wavelet transform modulus maxima 



(WTMM) at scale a is equivalent to the detection of edges of structures of size ~ a 6£ ] . 
This makes it a natural tool to search for the typical size at of mounds on the surface. 
Therefore, in the remainder of this chapter we use this particular wavelet. We investigate the 
average of the WTMM on the scale a 



n(a) := (AIr,\h](b,a))^ . (6.1) 

^ ^ \ v&MV , V6=WTMM ^ ^ 

This function has a pronounced maximum at a value which is the scale that contains 
the most relevant contributions to the surface morphology. On mounded surfaces, we expect 



Om oc Of. As an example, Q{a) of the surface of a simple cubic crystal is shown in figure 6.7. 



Figure S.6 shows the WTMM of the same surface at different scales including a = Om = 17. 
Clearly, at a = am the WTMM lines trace the contours of the mounds. On the contrary, 
at the small scale a = 5, the search for the WTMM yields structures at the flanks of the 
mounds which are small compared to the mounds themselves. The positions of the WTMM 
lines at a = 80 which is much greater than a^, are not related to the mounds at all. On other 
lattice structures and at different time, we obtain similar results. Consequently, a.^ is indeed 
a measure for the lateral size of the mounds. Since the convolution of the surface with the 
wavelet ^'o is the gradient of the surface smoothed with the gaussian fllter ^>o, ^l{am) is a 
measure for the height of the mounds. During the coarsening process, we have a^n ~ t^^^ and 
n{am) ~ t^- 
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a = ttr 



17 



a = 80 



Figure 6.6: WTMM at the surface of a simple cubic crystal {N = 512, 1^ = 15, Ik = 8, 
t = 3400 ML) at different scales a. The WTMM are shown in black. At a = am, the WTMM 
lines trace the contours of the mounds. At smaller a, the WTMM lines enclose details of the 
surface at the flanks of the mounds, whereas at large a the WTMM seem to be independent 
of the mounds. 
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Figure 6.7: Q(a) for the surface shown in 
figure 6.6. There is a pronounced maximum 
at a = Um = 17 where the WTMM lines en- 
close the mounds. 



This procedure has several advantages compared to the standard methods used in the 
literature: (1) The detection of the WTMM leads to an efficient suppression of noise and 
therefore eliminates the (noisy) intrinsic width. (2) Since only the most important lengthscale 
is considered, the results may not be corrupted by details of the surface morphology on small 
lengthscales. (3) Appropriate wavelets are well localized both in real space and in Fourier 
space. Therefore, this analysis combines the advantages of real space techniques, like the 
counting of mounds, and Fourier techniques, like the calculation of km- We have tested 
this method on various artificial test surfaces and found a clear superiority to the standard 
methods, especially in the presence of noise which should make it interesting for experimental 
investigations. 

6.5.2 Results 

In our simulations, the asymptotic scaling behaviour is obtained after a long initial transient 
of about 100 monolayers. Since after this transient only a comparatively small number of 
mounds is left on the surface, the measurement of exponents is complicated by finite size 
effects: due to the periodicity of the system, there is a self-interaction of the currents on the 
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mounds, if their size is no more small compared to the system size. We measure a ~ 1 for all 
Zfc > 1, a consequence of slope selection. At = 1 however, we measure smaller values which 
do depend on the crystal lattice: a = 0.85 on bcc and a = 0.94 on hep. A similar effect was 
observed in 101] on the simple cubic lattice. We remark that our results are independent 
of the method which was applied to measure the exponents. 
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Figure 6.8: /3, as obtained from the wavelet 
method, plotted as a function of the step edge 
diffusion length. +: bcc lattice, x: sc, □: hep, 
0: sh. Results from measurements of /3 from 
the surface width are identical within error 
bars which are on the order of magnitude of 
the symbol size. These have been estimated 
from the logarithmic fits; we expect the true 
errors to be larger due to systematical devi- 
ations. 



In figure 3. J, /3 is plotted as a function of Ik- Clearly, its value does not depend system- 
atically on the symmetry of the crystal surface. We find no significant deviations between /3 
on bcc and on hep (sc and sh). As already observed in ||l5| , |101| ], there is a dependence of (3 
on the step edge diffusion length. One obtains values ~ 1/4 at = 1 and higher values at 
greater Ik- The values obtained at large Ik are compatible with /3 = 1/3. As indicated above, 
we explain this transient as a competition between different coarsening mechanisms. At large 
Ik the merging of mounds should proceed mainly by the step edge diffusion driven transport 
of material into the gap, at small Ik it is dominated by noise assisted coarsening and/or bond 
energy driven coarsening. This competition might also explain why the exponents measured 
on bcc and hep are still slightly smaller than those obtained on the simple lattices even at 
values of Ik as large as 20, since the symmetry breaking of the mounds on these lattices halves 
the number of corners which are active in the coarsening process. 

All the simulation results reported above have been obtained in simulations with unhin- 
dered step edge diffusion. To further corroborate our ideas of the coarsening process, we have 
also performed simulations with a corner diffusion barrier. A particle may perform at most 
diffusion hops along a straight step edge, but diffusion around corners is forbidden. Then, 
we observe no symmetry breaking of the mounds, which now obtain regular polygonal shapes 
on all lattice structures. These are rotated by an angle of 45 degrees against the lattice direc- 
tions on cubic surfaces and 30 degrees on hexagonal ones. Due to the absence of aligned step 
edges, material transport by step edge diffusion is severely restricted under these conditions. 
Thus, we measure only small diffusion currents and observe no pronounced long-range order 
in the flux lines. Here, we find a ~ 1, /? ~ 1/4 on all lattices, independent of the step edge 
diffusion length. This slow coarsening in absence of the characteristic features of step edge 
diffusion driven coarsening — currents on mesoscopic lengthscales and symmetry breaking on 
bcc and hep — strongly supports the important role of this mechanism for fast coarsening with 
P = 1/3. 
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6.6 Conclusions 



In summary, we have presented a detailed investigation of the coarsening process in epitaxial 
growth. Our most important finding is that the crystal lattice has a considerable influence on 
the morphology of the growing surface which is by no means restricted to trivial symmetry 
effects. In spite of these differences, in the limit of large step edge diffusion lengths the 
mounds coarse according to power laws with universal exponents. We obtain /3 ~ 1/3 in the 
case of unhindered step edge diffusion. In the case of restricted step edge diffusion, /? 1/4 
is observed as it is for unhindered step edge diffusion with a small value of ^fc. 

Amar Q has obtained values of the coarsening exponents of bcc(OOl) surfaces which agree 
with our results within errorbars. He performed standard kinetical Monte Carlo simulations 
which consider the motion of many particles on the surface simultaneously. However, with this 
method only the simulation of the deposition of at most 300 ML was possible, and the initial 
formation of mounds took up to 100 monolayers. Our effective single particle dynamics, on 
the contrary, allows for the simulation of 2 • lO'^ ML such that coarsening can be investigated 
over two decades in time. Thus, our simulations give strong evidence that the mesured values 
of a and /? are indeed the dynamic exponents of the asymptotic regime. Conversely, the 
agreement of the results proves that the simplifications made in the effective single particle 
dynamics do not alter the values of the coarsening exponents. 

Our results contradict previous studies of the coarsening process using continuum equa- 



tions |6£, 110, 111] which predict slower coarsening on surfaces with a cubic symmetry. In 



these publications, equations have been studied which are invariant under the transformation 
h{x) — > —h{x). This symmetry reflects itself also in the up-down-symmetry of the solutions 
of the equations. Clearly, our simulation results break this symmetry and the arguments of 



1 95, lie, 111] cannot apply here. In figures 6^, QAh and 6^ the contours of the mounds 



and those of the valleys are clearly distinct. This is also true for the currents, especially the 
step edge diffusion current which dominates the coarsening behaviour in the limit of large Ik 
and unhindered step edge diffusion and determines the scaling exponents. We conclude, that 
the breaking of the up-down symmetry is a central feature of unstable epitaxial growth. Its 
precise effect on the coarsening dynamics clearly deserves further attention. 

However, we find it difficult to interpret our result in the context of an analogy to phase 
ordering processes, where the local slope of the surface corresponds to an order parameter. 
In this picture, the importance of the breaking of up-down symmetry leads to the conclusion, 
that now the stability of a domain wall between areas of equal slope is a complicated function 
of both the orientation of the wall and the order parameters on both sides of it. Additionally, 
in contrast to the simple lattices where there are only 4 (sc) respectively 6 (sh) stable slopes, 
on bcc and hep one would have to deal with a continuous order parameter. 

In any case, the behaviour of our models which implement the microscopic processes on the 
surface is governed by much more complex rules than those implemented in all the continuum 
models we are aware of. Differences are by no means restricted to some details, but have a 
fundamental influence on the behaviour of the system on large lengthscales. 
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Appendix A 

Transfer matrices and finite size 
scaling 



In principle, the transfer matrix method [|T^, 120 1 is a technique to calculate the free energy 



of one- dimensional lattice systems exactly in the thermodynamic limit. However, a system of 
arbitrary dimension can be treated formally as a one-dimensional system if it is decomposed 
in slices along an arbitrary direction. Each slice is then treated as one site in an infinite chain 
of slices. Since the transfer matrix method can consider only a finite, however large, number 
of possible states per site, the size of the slices must be finite. Therefore, in general we can 
investigate only semi-infinite systems which are finite in all directions but one. Only in some 
special cases, like the two-dimensional Ising model in the absence of a magnetic field, the 
limit of a real infinite system can be performed exactly. Otherwise, in two or more dimensions 
we must apply finite size scaling techniques to extrapolate to the thermodynamic limit from 
results obtained for semi-infinite systems. 



A.l The transfer matrix method 



Consider a chain of N identical slices of a system of arbitrary dimension. Let M denote the 
number of sites per slice. Each slice j can be in one of h different states Xj. We assume that 
there are interactions only between neighbouring slices. Interactions of longer but finite range 
can be considered by increasing the thickness of the slices. Then, the energy of the system 
can be written as 



7V-1 



H = H{Xn,Xi) + J2 H{Xj,Xj+i] 



(A.l) 



where we have assumed periodic boundary conditions. For simplicity, we assume without loss 
of generality that H{X,Y) = H{Y,X). It is always possible to rearrange the Hamiltonian of 
the system such that this property is fulfilled. The partition function is 



Af-l 



E 



exp 



(3H{Xn,Xi)\ J] exp -PH{X,,Xj+i) 



(A.2) 
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To calculate it, we map each of the n possible configurations of a slice to one of n orthonormal 
vectors, which span up a n-dimensional Hilbert space: 



1 if X = y 
else 



X < — >\X) such that {X\Y) = 

In the base of these vectors we define the h x h transfer matrix 

|X)exp[-/?^(X,y)] {Y 

\X),\Y) 



T 



(A.3) 



(A.4) 



which has the property that the trace of is equal to the partition function: 

TV[r^] = ^(Xi|r^|Xi) 

{Xi\T\X2){X2\T\X^)---{Xn\T\Xi) 



|^i),...,l^iv> 
Y exp 

|Xl>,...,l^iV> 

Z 



N-l 



-(iH{XN,Xi)\ J] exp \-l3H{X,,X,+i) 



Consequently, the partition function can be expressed in terms of the eigenvalues Aj of T. 
Since H{X,Y) = H{Y,X), T is symmetric and all Aj are real. We have 



Z = Tr [T^] =5^Af. 



(A.5) 



i=l 



This identity can be used to calculate the free energy per site / in the limit of infinite N. The 
theorem of Perron-Probenius states that the eigenvalue with the greatest absolute value is 
positive and nondegenerate. In the following, we will assume that the eigenvalues are sorted 
in descending order such that Ai > A2 > ... > Xa- Then, 



/ 



1 



InZ 



1 



In 



N- 



1 



InAi 



(A.6) 



Thus, we can calculate the free energy of a semi-infinite system by calculating the greatest 
eigenvalue of the n x n-dimensional transfer matrix T. Usually this is done numerically since 
n increases exponentially fast with the number of sites per slice. 

A. 1.1 Thermal averages 

The standard method in statistical physics to calculate the thermal average {q) of an arbitrary 
intensive quantity q uses a field h which couples to q. Then, 



(g) = lim ^ 
h^o dh 



(A.7) 



where fq is the free energy of the system with the modified Hamiltonian Hq = H + hNAiq. 
Of course, this method can be applied straightforwardly in the transfer matrix formalism by 
performing the differentiation with respect to h numerically. 
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There is however an alternative which can be apphed whenever q is weh defined as a 
function q{Xk) of the state of a single slice k. Then, since all slices are identical, the thermal 
average of q is the same in all slices and identical to the thermal average of q in the whole 
system. In a chain of N slices with periodic boundary conditions 

N-l 

= ^ E {Xi\T''-'\Xk)q{Xk){Xk\T''-''+'\Xi). (A.8) 

l^i>,l^fc) 

The transfer matrix T can be written as T = Y2i=i ' where jAj) is the normalized 

eigenvector to the eigenvalue Aj. Further, we introduce the matrix 

q = l^> li^) (^1 (A.9) 

\x) 

which will be denoted as the operator of q in the following. Then, 

where we have used that (AjjAj) = 5ij and l^i) — 1- Thus, the thermal average 

of q can be calculated easily from the normalized eigenvector |Ai). 

In particular, the probability P{X) to find a slice in state X is the thermal average of 
nx which is defined as one, if the slice is in state X and zero otherwise. Its operator is 
hx = \X) {X\. Using equation [A.10| , we find 

P{X) = (Ai|X) {X\Xi) = \{Xi\X)f . (A.ll) 

The square of every component of the normalized eigenvector to the greatest eigenvalue of 
T, represented in the base of the vectors \X), is the probability to find a slice in a particular 
state. 



A. 1.2 The correlation length 



If we measure an arbitrary quantity q simultaneously in two different slices, due to the inter- 
actions between the slices these values are not statistically independent. This fact is measured 
by the correlation function 



g,{\l-k\):= {q{Xk)q{Xi)) - {qy 



(A.12) 



Nonzero values indicate correlations between q{Xk) and q{Xi). The thermal average of 
q{Xk)qiyXi) in a chain of N slices can be expressed in terms of the transfer matrix: 



N-l 



{q{Xk)q{Xi)) 



- q{Xk)q{Xi)exv[-(5H{XN,Xi)\\\ex^\-PH{Xi,Xi+i 
^ {Xi\T^-^qT^-^qT^-^+^ \Xi) . (A.13) 



1^1 ) 
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For simplicity, we have assumed that I > k. We continue by expressing the transfer matrix in 
terms of its eigenvectors and eigenvalues and performing the limit N ^ oo. We find 



{q{Xk)qiXi)) 



^j2xt'{X^\ci\X,)X'r^{X,\q\X^)Xr^' 



{Xi\q\X^f + Y, 



l-k 



l(Ai|g|A,)r 



(A.14) 



The first term in equation A.14 yields (q) . If / — /c ^ 1, the sum is dominated by the term 
where the absolute value of Xj is maximal and the matrix element {Xi\q\Xj) 7^ 0. In the 
following, we denote this eigenvalue as Am and the corresponding eigenvector as |Am,)- This 
yields the following asymptotics for the correlation function: 



g,i\l-k\)^\{X, 



|A™)|M'^ 



\l-k\ 



(A.15) 



If Xm is positive, gq{\l — k\) decays exponentially. On the other hand, if Am is negative, 
gq{\l — k\) oscillates with an exponentially decaying amplitude between positive and negative 
values. This indicates that the values of q in neighbouring slices are anticorrelated. In both 
cases, the amplitude decays on a typical lengthscale 



In 



Ai 



(A.16) 



which is called correlation length. 



A. 1.3 Symmetries of the Hamiltonian and vanishing matrix elements 

The transfer matrix formalism is based on the assumption that the Hamiltonian of the inves- 
tigated system is invariant under translations of the system in the infinite direction. However, 
frequently both the Hamiltonian and the operator of a quantity q fulfil additional symmetry 
properties. These can be used to predict which of the matrix elements (Aj| ^|Aj) are zero. 
Let there be an operator U which maps each configuration of a slice uniquely to another one. 
We assume that the energy of the system is invariant under this transformation. Since the 
transfer matrix has the same symmetry properties as the Hamiltonian, we have 

U^U=1 (A.17) 



UTU^ = U^TU = T. 



(A.18) 



Further, we assume that the operator q is either symmetric or antisymmetric under the trans- 
formation U such that 

UqU^ = Xqq where Xq e {-1, 1}. (A.19) 



Applying equations A.17 and A.l£ to the product of the transfer marix with its eigenvector 
I A, ) we find 



UTlXi 



UTU^UlXi) = TUlXi) = X,U \Xi 
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(b) 




Figure A.l: Sketch of a square lattice system as discussed in section |A.1.4 . Panel (a) shows 
the decomposition of the system in lines with L sites each. Panel (b) shows the unit cell of 
the lattice. Sites are shown as black circles. Possible interactions are symbolized by lines. 
Note, that each of the sites belongs to four unit cells. The nearest neighbour bonds belong to 
two unit cells. Therefore, in the calculation of h{v,w,x,y), the energies of these objects have 
to be weighted with the factors shown in the figure to avoid multiple counting of energies in 



equation A. 24 



Consequently, U W) is an eigenvector to the eigenvalue of the transfer matrix. In case of 
a non-degenerate eigenvalue Aj, this condition implies 



U \Xi) = Xi where Xi G {-1, 1}, 



(A.20) 



since due to equation [A.17| the multiplication with U leaves the norm of |Aj) invariant. If 
Xi is degenerate, the eigenvectors of T are not uniquely determined. However, it is always 



possible to construct a set of orthonormal eigenvectors such that equation A.20 is fulfilled. 
Vectors with Xj = 1 are symmetric under the transformation U. Xi = — 1 characterizes an 



antisymmetric vector. Using equations A. 17, A. 19 and A.20, we find 



{Xi\q\Xj) = {X^\U^UqU^U\XJ) =x^XjX4{>^^\Q\>^J) ■ 



(A.21) 



For a nonzero matrix element (Aj| q\Xj), this equation can be fulfilled only if XiXjXq = 1- 
Consequently, the matrix elements of a symmetric operator q, where Xg = 1 vanish whenever 
the eigenvectors |Aj) and \Xj) have different symmetry properties under the transformation U. 
Conversely, for an antisymmetric operator (xq = ~1) matrix elements between eigenvectors 
of equal symmetry must be zero. 



A. 1.4 Application to two-dimensional lattice systems 

Now we apply this formalism to a special class of models which are defined on a two- 
dimensional square lattice. Each lattice site can be in one of n different states. There are 
interactions only between nearest neighbour and diagonal neighbour sites. The models dis- 
cussed in section belong to this class. We consider a rectangular system with an extension 
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of L lattice constants in the x-direction and lattice constants in the y-direction as shown 
in figure [A.l| a. The boundary conditions are assumed to be periodic. 

To investigate the limit — > oo of an infinite strip with the transfer matrix method, we 
decompose this system into slices which consist of a single line of sites in the x-direction. 
The configuration X of such a line is determined by the states of the L sites (xi,X2, ...^xi). 
Since there are n possibilities for each of the Xj, we have n = different configurations of 
a line. The transfer matrix formalism requires a mapping of these configurations to a set of 
vectors such that {X\Y) = 5xx,yi ■ ^X2,y2 ' ' ' ^XL,yL any two vectors \X), \Y) which represent 
configurations X, Y. For this purpose, we introduce real n^-dimensional unit vectors 

\^) = (^^1,(;{X),^2,((X), ■■■,SnL,({X)) J (A.22) 

where C(^) is a function which maps each configuration uniquely to an integer € {1, ...,n^}. 
To define such a function, we assume without loss of generality that the Xj can have values 
from to n — 1. Then, we can read each configuration of a line as a number G {0, ...,n^ — 1} 
written in a positional notation with radix n. Consequently, 

L 

C(A) = l + ^x,n^-i (A.23) 
1=1 

has the desired properties. It is convenient to use the notation \X) = |xi, ...,xl) which shows 
the configuration of the line that is represented by \X) directly. Note, that the Xj are not the 
components of \X)\ 



The Hamiltonian of the system can be written in the form of equation |A.l| , where 

H{X,Y) = h{xL,xi,yL,yi) + ^ /i(xi, Xj+i, y^, yi+i). (A.24) 

1=1 

Here, h{v, w, x, y) is the total energy of a unit cell as shown in figure |A.l| b. Inserting equation 
A.24| into equation |A.4| , we obtain the transfer matrix 



L-l 

^= X] \^^'-^^L)t{xL,Xi,yL,yi)Y[iiXi^^i+l^yi^yi+l){yi^-^yL\, (A.25) 
\X),\Y) i=l 

where we have introduced t{v,w,x,y) = exp{—f3h{v,'w,x,y)). 

We are interested in the extremal eigenvalues of T and the corresponding eigenvectors. 
Since T is a x n'^ matrix, the number of its elements, n^^, increases extremely fast with L. 
Therefore, any approach where the transfer matrix is stored in computer memory is limited 
to very small strip widths L. However, there are algorithms for the computation of extremal 
eigenvalues which do not need the matrix itself but only a method to compute the product 
of the investigated matrix with an arbitrary column vector. One of these is the Lanczos 



algorithm which we use in this work. We refer the reader to the numerical literature |114| for 



a description of this method. In this context, it is important to note that the transfer matrix 



can be written as a product of L matrices which are extremely sparse [78, 12C]: 

T = TlTl-1 • • • Ti (A.26) 
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The n^'^'^ x matrix 

7i = XI X] \yi^^i^-^^L,yL+2)t{yL+2,yi,XL,xi){xi,...,XL\ (A.27) 
\x) yi,yL+2 

maps the n^-dimensional space which is spanned up by the vectors {\xi, ...,xl)} to a n^"*"^- 
dimensional space. An orthonormal base which spans up this space can be constructed using 
a mapping of the configurations of a line of L + 2 sites to vectors similar to that defined in 



equations A. 22 and A.23 . The matrix elements are then defined in terms of the configuration 
both of the line of L sites and the line of L + 2 sites. Note, that Ti has only • nonzero 
elements. Similarly, we have L — 2 matrices of dimension n^^^ x n^"^^ 

Ti = 1^1' ■■■^Xi-i,yi,Xi+i, ...,XL+2) t{xi-i,yi,Xi,Xi+i) {xi, ...,xl+2\ (A.28) 

\X) Vi 

where i = 2,...,L — 1. Each of these matrices has n ■ n^^^ nonzero elements. Finally, we have 
the X n^"*"^ matrix 

Jl = X] 1^1' ■■■,XL-l,XL+2) t{xL-l,XL+2,X2,XL+l) {xi, ...,Xl+2\ (A.29) 

which has nonzero elements. It is easy but lengthy to verify that the product of 

these matrices yields the transfer matrix. Using a sparse representation, we need to keep 
only the nonzero elements to store the Tj in computer memory. The Tj altogether have 
^2 -\- (^L — 2) n] nonzero elements, which is significantly less than the n^^ elements of T 
itself. Since the product of the transfer marix with an arbitrary vector \v) can be performed 
by successively multiplying \v) with the Ti, we can save a great amount of memory and CPU 
time by using this factorization. On our computer equipment, we can deal with strip widths 
up to L = 18 for n = 2 and L = 10 for n = 3. 

A. 2 Finite size scaling 

Using the transfer matrix formalism, the thermodynamic properties of semi-infinite systems 
can be calculated exactly. However, we are interested in the thermodynamic limit of a system 
which is infinite in all directions. In particular, we want to study the phase transitions 
of the infinite system. A phase transition is an abrupt change of the order of the system 
in dependence of a parameter like temperature. At continuous transitions and first order 
transitions which are investigated in this thesis, the correlation length diverges to infinity. 



Additionally, there is a singularity in the free energy. From equation [A.lq , it follows that a 
diverging correlation length requires a degenerate greatest eigenvalue of the transfer matrix. 
Consequently, there are no phase transitions in semi-infinite systems since there the greatest 
eigenvalue of the transfer matrix is always nondegenerate. On the other hand, it is reasonable 
to assume that we can draw conclusions on the infinite system from results obtained for a 
semi- infinite system with sufficiently large slices. The theory of finite size scaling describes 
how this can be done. The procedure is different for continuous phase transitions and first 
order transitions. 
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A. 2.1 Continuous phase transitions 

In continuous phase transitions, the order parameter m which measures the order of the system 
decreases continuously to zero as the critical temperature Tc is approached from below. There 
is no latent heat and no interface tension between the low temperature phase and the high 
temperature phase. The description of this kind of transition is based on the observation that 
the system is self-similar at the critical point, i.e. the statistical features are invariant under 
rescaling. The method which is presented in this section has been suggested by Nightingale 
fn\ , [78| as "phenomenological renormalization group theory" . 

We introduce the parameter e := T — T^. Additionally, we consider an ordering field h 
which couples to the order parameter m. Then, in the infinite system the critical point is 
at e = /i = 0. For simplicity, we assume an equal system size L in all finite directions. We 
divide the system into hypercubes with a linear dimension of b lattice sites. For small enough 
e, h we have 6 <C ^. We map each block to a single site in a system of size L/b which might 
be done e.g. by a majority decision among the sites in the block. Following renormalization 
group theory, the free energy of the interactions between the blocks can be approximated by 
the free energy of the rescaled system at a different temperature Cb and field hh- Then, the 
total free energy per block is 

6'^/(e, h, 1/L) = b'^Qbie, h, 1/L) + /(e^,, h, b/L). (A.30) 

d is the dimension of the system. The first term on the right hand side of this equation is the 
free energy of the interactions between the sites in the block. Since b is finite, it is analytical 
even at the critical point where e = h = Q and L = oo. All free energies have been expressed in 
terms of free energy densities per site. It is convenient to define the free energy as a function 
of the inverse system size. Then, at the critical point all arguments of / are zero. Clearly, 
the correlation length in the rescaled system is by a factor l/b smaller than in the original 
system such that 

^{e,h,l/L) =bC{eb,h,b/L). (A.31) 

A simple consideration yields the functional dependence of ef, and hj, on e and h. We 
iterate the rescaling procedure by dividing the rescaled system into blocks of size b'. Clearly, 
this should yield the same result as a division of the original system into blocks of size bb': 
e6'(ef)(e)) = ^b'bi^) and hy{hb{h)) = h^'hih). This group property is fulfilled if 

ef,(e) = by^e 

h{h) = b^^h ^^■'^^> 

with scaling exponents yx and yh- These exponents characterize the universality class of the 
model. Note, that formally the inverse system size enters the free energy and the correlation 
length like a field which drives the system away from the phase transition. The corresponding 
scaling exponent is one. 

One can show that the free energy can be decomposed into a singular and a regular 
part / = /sing + /rcg- /rcg IS analytical even at the critical point. Close to the phase transition, 
the behaviour of the system is dominated by the singular part /sing which fulfils 

&^/smg(e, h, 1/L) = fsingiby^e, b^^h, b/L). (A.33) 



^In principle, the temperature might be replaced by any other parameter which can be tuned to reach the 
critical point, e.g. a chemical potential. 
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The second derivatives of this equation with respect to e and h, respectively, yield similar 
scaling forms for the specific heat cl and the susceptibility xl of m. From these, we obtain 
how the singularities oi cl and xl build up as L goes to infinity: 



XL oc hh{o,o,i/L)^L^y'^-^. 



(A.34) 
(A.35) 



Using equations A. 31 and A. 32, we obtain the scaling behaviour of the inverse correlation 
length K := 

bK{e, h, 1/L) = K{by^e, b^^h, b/L). (A.36) 

Consequently, at e = /i = the inverse correlation length, its derivative with respect to 
temperature and its second derivative with respect to the field h behave like 

-1 



«;(0,0, 1/L) 
«.(0,0,1/L) 
Khh{0,0,l/L) 



(A.37) 
(A.38) 
(A.39) 



In equation A.39, we have assumed h ^ —h symmetry of the investigated model such that the 
first derivative of k with respect to h vanishes. These properties are sufficient to determine 
the critical temperature and the scaling exponents ut and yh- To this end, we measure the 
inverse correlation lengths at zero field for two different strip widths L^^^ and L^^^ 

as a function of temperature. Due to equation A.37 , Tc is the temperature where k^^^L^^^ = 
^{2)^(2)_ Q22ce Tc is known, yx and y^ can be determined from equations A.38 and A.39| . At 



Tc, we have 



In 



In 



L{2) 



and Vh 



2 U'^'41 



In 



L(2) 



(A.40) 



A. 2. 2 First order phase transitions 



At a first order phase transition, there is a discontinuity both in energy and the order pa- 
rameter. The energy difference between the high temperature phase and the low temperature 
phase is called latent heat. Additionally, there is an interface tension between the phases. 
Since the surface of a phase and its volume scale differently if the rescaling procedure intro- 



duced in chapter A. 2.1 is applied, this interface energy forbids self-similar properties of the 
system. It is connected with the phenomenon of metastability: The low (high) temperature 
phase is insensitive to small perturbations even at temperatures above (below) the transition 
temperature. Due to the interface tension, small droplets of the stable phase in the unstable 
phase are unfavourable. At the transition temperature, the system will spontaneously order 
in one of the phases unless this is impossible due to some constraint, e.g. a fixed particle 
number. In that case, the coexisting phases are separated and exist in different parts of the 
system. 

The description of these phenomena requires an extension of the concept of free energy. 
Conventionally, the free energy is defined as a thermodynamic potential which characterizes 



the state of the system. In the theory of phase equilibria [92|, the concept of the free energy 
of a phase is used. In the thermodynamic limit, the free energy of the system is then assumed 
to be the sum of the free energies of all coexisting phases. In the following paragraph, we 
discuss the connection of this idea with the transfer matrix formalism. 
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Phases and asymptotically degenerate eigenvalues 



Due to the dependence of the correlation length on the eigenvalues of T in a finite system 
(equation [A.(]| ), one might expect that in the thermodynamic limit the greatest eigenvalue of 
the transfer matrix becomes degenerate at a phase transition. However, the transfer matrix 
of an infinite system is ill-defined. Since energy and free energy are extensive, H{X, Y) and 
the free energy per slice grow proportional to the number of sites per slice, M. Therefore, 
some elements and eigenvalues of T go to infinity as L — > oo. To avoid this mathematical 
problem, we define the quantities 



fi-.-- 



-1 



/3M 



In I A,- 1 for 



1, ...,n. 



(A.41) 



The eigenvalues Aj are assumed to be sorted in descending order. For L ^ oo, fi converges to 
the free energy per site of the infinite system. Consider a system at a phase transition. Quite 
generally, we assume that there are z/ positive eigenvalues Ai > A2 > ... > Ai/ and /.t negative 



eigenvalues An < An_i < ... < Xn- 



-i-i+i 



such that 



Ai 

A,: 



-PM (/i 



/.) 



(A.42) 



for i < 1/ and i > h — ^. In the following, we will denote these eigenvalues as "asymptotically 
degenerate". Each of the might be the diverging correlation length of an appropriately 
chosen quantity. 

We investigate the influence of a perturbing field h which couples to the extensive quantity 
NAiq. q is an arbitrary intensive quantity the operator of which is q. This field drives the 
system away from the phase transition at /i = 0. In the remainder of this paragraph it is 
assumed that \h\ is so small that terms of order /i^ and higher orders in h can be neglected. 
Variables which refer to the perturbed system are marked with a prime. The transfer matrix 
is 



T' = Yl i^)«^p 
= E 

\X),\Y) 

(3hM 



p{H{X,Y) 



hM 



{q{X)+q{Y)) 



exp 



-PH{X,Y) {Y\ 



T + 



(qT + Tq). 



(A.43) 



According to perturbation theory, for small L the greatest eigenvalue of T' is A'^^ = (Ai| T' \ Xi) + 
0{h'^). However, for large L the gap between the greatest eigenvalue of T and the asymptot- 
ically degenerate eigenvalues becomes smaller. For any small but finite h, there will be some 
value of L where the perturbation induced by h is comparable to the separation between the 
eigenvalues. In this case, the eigenvector {X^) to the greatest eigenvalue of T' will be the lin- 
ear combination of the eigenvectors to the asymptotically degenerate eigenvalues of T which 
maximizes (A'^| T' |A'^). We make the ansatz 



lA'i) 



U+fJ. 

E 

k=l 



Ck 



\K{k)) 



(A.44) 
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where Ylkt^i ^1 = 1 and 



zik) 



k if k < v 

k + h — V — II else 



(A.45) 



enumerates the asymptotically degenerate eigenvalues. Using equations |A.44| and |A.43| , we 
obtain 



(A;|r'|A'i> = Ai ckci 

kl=l 



A 



Ai 



4,/ + 



(3hM 



Ai 



Ai 



Due to equation A.42 , we have limi^ooAj/Ai = sign(Aj) for all asymptotically degenerate 
eigenvalues. Using this fact and equation |A.6| , we obtain after some lengthy algebra that in 
the limit L — > oo 

W _ / /a = /l - %max if /l > ( K Ap.\ 



where gmax (gmm. 



is the greatest (smallest) eigenvalue of the matrix 



Q 



( 

V 



Q 



n—fi+l,n—fj,+l 



Qn 



Qn,n 



h—fi+l 



n—fi+l,n 



Qn,n / 



(A.47) 



Here, we have introduced Qij = lim^^^oo (Aj| q\^j)- This result is reminiscent of degenerate 
perturbation theory. Note, that there is no mixing between positive and negative eigenvalues. 



Using equation |A.7| , we find that there might be different expectation values of q at the phase 
transition depending on whether the limit /i — > is performed for a positive or a negative 
field: 



{qa 



hm — 

h^o+ oh 



lim 



dh 



(A.48) 
(A.49) 



Clearly, this is what one would expect at a first order phase transition where there is a 
discontinuity in q. Obviously, this is the case if q is an order parameter. On the contrary, at 
a continuous phase transition we have ^max = 'Zmin for an arbitrary quantity q. This implies 



that Q has only one v + ^ -fold degenerate eigenvalue. Figure A. 2 shows a sketch of the free 
energy /{. As expected, at a first order phase transition there is a cusp in the free energy. At 
a continuous transition, the first derivative of /{ with respect to h is continuous. 

Consider a system with a first order phase transition and z/ + /i = 2 asymptotically de- 
generate eigenvalues. In the absence of constraints, at /i = the symmetry of the system is 
spontaneously broken. The system is either in a phase where {q) = Qmax (denoted as phase 
(a) in the following) or a phase where {q) = gmin (termed phase (b)). For /i > 0, phase (a) is 
thermodynamically stable, for /i < (b) is the stable phase. If the concept of the free energy 
of a phase can be applied, fa should be the free energy of phase (a). Conversely, the free 
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Figure A. 2: Sketch of the behaviour of /{ at a phase transition where + /x = 2 in the 
hmit L ^ oo. Panel (a): First order phase transition. The sohd hne shows the free energy of 
the system. At the phase transition (/i = 0), there is a discontinuity in df[/dh. The dashed 
hues show the continuation of the free energies of the phases into the regions of metastabihty. 
Panel (b): Continuous phase transition. At /i = 0, the first derivative of /{ is continuous. 



energy of phase (b) should be /;,. We expect that both phases are metastable in the region 
where they are thermodynamically unstable. Since there is no singularity in the free energy 
of a phase at the transition (in contrast to the free energy of the system) , the free energies of 
the phases in the metastable region should be given by the continuation of fa and /;,. 

Generalizing this idea, we find that in the vicinity of a first order phase transition the 
free energies of the metastable phases are functions of those eigenvalues of the transfer matrix 
which are asymptotically degenerate at the phase transition. The number of these eigenvalues 
equals the number of coexisting phases at the transition itself. For example, at a triple point 
we expect u + /i = 3. 



Application to semi-infinite systems 

At continuous transitions, the inverse system size enters the free energy like a field which drives 
the system away from the phase transition. According to Privman, Fisher and Schulman 
|]89| , |90t| , a similar analogy holds in the vicinity of a first order transition. They argue, that 
the fi of a semi-infinite system can be calculated from measurends of an infinite system via a 
perturbation theory similar to that presented in the preceeding paragraph. The finiteness of 
the system size plays the role of a perturbing field. Consequently, the free energy per site in a 
semi- infinite system is the smallest eigenvalue of a {v + fi) x {f + fi)- dimensional free energy 
matrix which has the same structure as Q in equation A. 47] . Its diagonal elements are the 



free energies per site of the phases which coexist at the transition. The off-diagonal elements 
determine the finite size corrections. Since there is no mixing between positive and negative 
eigenvalues, the free energy matrix can be diagonalized blockwise. 

In the models which are discussed in this thesis, we have to deal with lines of triple points, 
lines of points where four phases coexist and, in the model of the CdTe(OOl) surface with Te 
dimerization, a point where five phases coexist. In our investigations, the system is arranged 
such that the repulsive interactions are in the infinite direction. Then, neighbouring slices 
are anticorrelated which implies negative asymptotically degenerate eigenvalues. We find, 
that whenever the number of coexisting phases is < 4, there are at most two positive and 
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h 



Figure A. 3: Panel (a): Schematical behaviour of the fi at a finite strip width L. The 
dashed lines show the free energies fa, fb of the phases of the infinite system. Panel (b): 
Sketch of a semi-infinite system at the phase transition of the infinite system. White and grey 
areas symbolize two coexisting phases. 



two negative eigenvalues which are asymptotically degenerate. In the following, we restrict 
ourselves to this case. The point where 5 phases coexist can be found from the intersection of 
lines of other phase transitions. Due to the theorem of Perron-Frobenius, /i is the smallest 
eigenvalue of the positive block of the free energy matrix which has the form 



fa A 

A fb 



(A.50) 



Its eigenvalues are 



/i,2 = -^111^1,2 = 2 (/- + f") T 2 V (-^^ - f''^' + (^-51) 

At the phase transition, we have fa = fb- Then, the inverse correlation length is = 
In IA1/A2I = 2/3A. Further away from the phase transition, where \fa-~ fb\ ^ A, fi approaches 
to the minimum of and (figure |A.3| a) . 

Since A is a function of the correlation length at the phase transition, we can determine 
the functional form of the dependence of the finite size corrections on L from an estimation 
of ^. In a semi-infinite system, there is no spontaneous symmetry breaking at the transition 
temperature of the infinite system. Instead, the system consists of alternating domains of the 



coexisting phases. A sketch of such a system is shown in figure A.Sb. The correlation length 
is proportional to the typical distance between two domain walls. Due to the interface tension 
a, the formation of such a domain wall costs a free energy aLa, where a is a dimensionless 
number 0{1). Consequently, the probability for the formation of a domain wall is proportional 
to exp{—(3aLa). Therefore, we have 

^ oc exp {(3aLa) =^ A oc exp {—AL) where A = f3aa. (A. 52) 

Since the behaviour of all other quantities depends on the free energy, an exponential decay of 
the finite size corrections should be a generic feature of first order phase transitions. However, 
frequently one encounters the phenomenon of a tricritical point at a temperature Ttri. At 
Ttri, the nature of the phase transition changes from first order at T < Ttri to a continuous 
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transition at T > Tti-i. Conseqently, at Ttri we expect a crossover from exponential finite size 



scaling to the linear scaling of the continuous transition 1 93 ] . Unfortunately, at feasible values 
of L, this crossover is smeared out over a finite temperature range which complicates finite 
size scaling in the vicinity of the tricritical point. 

To determine the locus of the phase transition, we consider an intensive quantity q and 
a field h which couples to NAiq. We assume the phase transition to be at /i = 0. Since 
equation 



A. 51 holds also in the presence of a field, we have 
dfi,2 



dh 



dfa Oft 
dh dh^ 



dfa_dfb 

dh dh 



fa fb 



(/a-A)'+4A2 



At the transtion, we have fa = fb- Consequently, 

dfi,2 _ 



(9)1, 



2 •- 



1™ au 

h.-»o dh 



(dfa 




dfb 


V dh 


h=0 


dh 



h=0 



(A.53) 



Using equation A.10| , we can calculate (q)-^ 2 from the eigenvectors |Ai) and IA2) of the transfer 
matrix. Then from equation A.53 at the phase transition we have 



(All q |Ai) = (A2I g|A2) . 



(A.54) 



Note, that this identity is independent of the order of the phase transition and the functional 

form of the finite size corrections. 

Once the locus of the transition is known, the values of q in the phases can be determined 

.(/) 



from the matrix elements {Xi\q\Xj) =: Ql j ■ Consider the matrix 

/ Q'/J ■ 



)(/) 



Ql'l 



\ 



^ h,n—fi+l 



Q 



if) 



(A.55) 



Since lim^^oo Q^^"^ = Q-, the eigenvalues q^^"* of Q^f^ converge to the eigenvalues qi of Q in 
the limit L — > 00. Due to equation A. 52, we expect 



q\^^ =qi + Be^Y>{-AL) 



(A.56) 



Therefore, qi and the parameters B can be determined from transfer matrix calculations 
with three different strip widths. 

Methods similar to that presented in this section have been suggested in the literature for 



the investigation of the spontaneous magnetization |41] and latent heat [42| of the Ising model 
and for the calculation of coverage discontinuities in interacting hard square lattice gases 1 14 ] . 
However, to our knowledge a consistent presentation of this subject in the general case has 
not been published yet. 



103 



Appendix B 

Continuous time Monte Carlo 
simulations 



A large number of Monte Carlo algorithms is based on the mathematical concept of Markov 
chains [^]. In all the simulations presented in this thesis, such methods were used. Consider 
a system with a finite set S of possible states s. A Markov chain is a stochastic process 
in discrete timesteps which is characterized completely by a set of transition probabilities 
{W<j^<j/|s, s' S S}, where Ws^s' is the conditional probability to find the system in state s' at 
step n + 1 given that it is in state s at step n. 

In equilibrium simulations one is interested in the calculation of averages over an ensemble 
of states distributed according to some probability distribution function -P(s). Of course, the 
most important example is the Gibbs distribution 

exp[-H{s)/ikT)] 
^ ' Es^^P[-H{s)/{kT)r 

where H{s) is the energy of the system in state s, k is Boltzmann's constant and T is tem- 
perature. The idea is to construct a Markov chain, where time averages over this sequence of 
states converge to ensemble averages over P{s) in the limit of infinite time. The first condition 
for this to be the case is that every state s of the system with P{s) > can be reached within 
a finite number of steps from every initial configuration. A dynamics with this property is 
called ergodic. The second condition is that P{s) is stationary under the dynamics. Consider 
an ensemble of Markov chains. Let Pn{s) denote the probability of a randomly chosen member 
of such an ensemble to be in state s at time step n. The change of this probability in the next 
step is then given by 

P„+l(s) - Pn{s) = Ws'^sPnis') - Ws^s'Pnis). (B.l) 

s'^s s'^s 

The first sum is the total probability to enter state s in step n, the second sum is the total 
probability to leave state s. A distribution is stationary, if the left hand side of this equation 
is zero. In particular, this is the case if the transition probabilities fulfil a detailed balance 
condition 

Ws^s'P{s) = Ws'^sP{s'). (B.2) 
In case of a Gibbs distribution, this condition yields 

^-'-expf-^ifl^V (B.3) 



VF,'^. V kT 
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Since there is no unique solution of equation B.3, there is some freedom in the choice of 
transition probabiUties in simulations of thermal equilibrium. In particular, there is no need 
to choose the Wg^g' in agreement with the physical dynamics of the system. 

In the simulation of systems far from equilibrium like crystal growth or sublimation one 
is interested in the evolution of the state of the system with time. If we look at the state of 
the system on a coarse-grained scale in discrete time intervals of appropriate length, we can 
describe its evolution by a Markov process. The transition probabilities depend on the chaotic 
microscopic motion of atoms (section [1.41 ). Therefore, the Wg^s' must reflect the physics of 
the investigated system. However, in the absence of external forces any physical system will 
finally reach thermal equilibrium. Therefore, any reasonable choice of the dynamics of the 
isolated system must fulfil at least the conditions of ergodicity and the stationarity of the 
Gibbs distribution. The simplest way to ensure the latter condition is to make an ansatz for 
the Wg^gi in agreement with the detailed balance condition B.3 . 

It is favourable to express the transition probabilities in terms of events which modify the 
system. Such events might be e.g. diffusion hops of atoms on a crystal surface or spin flips in 
an Ising model. In general, let there be a system where M events i = 1, ...,M are possible. 
There is a rate Pi{s) for every event, which is a function of the actual state s of the system. 
We introduce the transition probabilities 



Wg 



•Af pi{s) 



1 JTpo 



if s 7^ s' and i{s 
if s = s' 
otherwise 



s ) exists 



(B.4) 



where po is a constant such that po ^ max,, pi. Event z(s — > s') changes the state of the 
system from s to s'. These Wg^s' fulfil the conditions Wg^s' ^ and Yls' ^s^s' = 1 for 
probabilities. Detailed balance is fulfilled if Pi(^g_^gi^/ Pi(gi^g-) = P{s) / P{s'). 

These transition probabilities describe the dynamics of a stochastic system with indepen- 
dent events if the probability that two or more events occur simultaneously in the physical 
time interval At between subsequent Markov steps can be neglected. The probability of a 
random event with rate Pi{s) to happen in a time interval At is pi{s)At. Identifying this with 
the corresponding transition probability we find 



At 



1 



Mpo' 



(B.5) 



Clearly, the probability for two or more events in At decreases ~ 1/Po Po increased. In 
continuous time algorithms, we can perform the limit pQ ^ oo exactly, such that there are no 
problems due to simultaneous events. 



B.l Fundamental algorithms 



In the Markov chain defined in equation B.4, there are two different classes of steps: On 



the one hand, there are steps in which the state of the system changes. In the following, 
we will denote these as "modifying steps". On the other hand, there are steps in which the 
state of the system remains the same (non- modifying steps). Especially at low temperature 
the probability Wg^g of a non-modifying step usually is high. This leads to a considerable 
slowing down of simple techniques like the Metropolis algorithm which simulate the Markov 
chain directly. 
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Continuous time algorithms [|T^, ^ avoid this problem. The basic idea is to jump over 
the non-modifying steps. In every iteration of the algorithm one event is selected randomly 
with the conditional probability Pi{s) that event i is performed in a Markov step given that 
this is a modifying step. Then, the number Up of Markov steps which pass before the next 
modifying step is drawn randomly according to the probability of non- modifying steps. This 
procedure yields the correct statistics of the Markov chain, if quantities measured in the actual 
state of the system are weighted with Up or, equivalently, with the corresponding physical time 
interval. According to equation B.4, the probability of a modifying step is 

M 



Pm{s) = l-Ws^s = ^ where = V (B.6) 



i=l 



Using this equation, we can calculate Pi{s). The probability that event i is performed in one 
Markov step is 

Npo R[s) 

Consequently, 

The probability of a specific value of Up is the probability of a sequence of np — 1 non-modifying 
steps followed by one modifying step which is 

P{np, s) = {l- PMr^-'P^is). (B.8) 

To perform the limit po — > oo, we express Up in terms of the physical time interval r = 
UpAt = np/{J\fpo) which passes between subsequent modifying steps. Then, the probability 
distribution function of r is 

/ E(s)\-^^°^ 

P{t, s)dT = lim 1 - -f^ R{s)dT = R{s) exp {-R{s)t) dr, (B.9) 

po->oo Y J\l Po J 



where we have used equations B.6 and 

Algorithm B.l Continuous time Monte Carlo simulation 

1. Draw an event i randomly with probability Pi{s) and perform it. This changes the state 
s of the system. 

2. Calculate the Pi{s) and R{s) in the new state of the system. If the event changes the 
state of the system only locally and the range of interactions is short, pi{s) changes only 
for a few i . 

3. Draw a uniformly distributed random number x G ]0, 1] and increase the system time by 
T = — lnx/i?(s). T is a random number distributed according to equation B.§. In the 
calculation of thermal averages, weight the quantities measured in the current state of 
the system with r. 

4. Co to step 1 unless some condition to terminate is fulfilled (e.g. a maximum number of 
steps is done). 
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Figure B.l: Sketch of a complete binary tree. 



To apply this algorithm, we need fast methods to perform the random selection of events in 
step 1 and the update of the rates in step 2. We first introduce algorithms which perform 



both tasks in C'(log2AA) time [17j. They can be applied in the general case where there is a 
unique rate for every event. 

We use a complete binary tree [^3| to store the list of the rates of all events. Formally, a 
binary tree is defined as a finite set of nodes that either is empty, or consists of a root node 
and the elements of two disjoint binary trees which are called the left and the right subtree of 
the root. The root nodes of the subtrees of a node are commonly denoted as its "children". 
Conversely, the node a node is child of is called its "parent" . 

Usually, a binary tree is drawn with the root node on top. Lines connect every node with 



its children. Such a sketch of a binary tree is shown in figure B.l. There is a simple way to 
enumerate the nodes of a binary tree. The root node gets index 1. In general, node k has the 
left child 2k and the right child 2k + 1. Then, the parent of node A; is [A;/2j. A binary tree 
is called complete if this enumeration maps its nodes to a set of consecutive integers. The 



tree shown in figure B.l is such a complete binary tree. We use this enumeration to store a 
complete binary tree in a linear array in computer memory. 

In our Monte Carlo simulations, we use a complete binary tree with N nodes, where M 
is the maximum number of events. Each node k is attributed to a specific event. There are 
three real numbers for each node: Wk = Pk{s) is the rate of event k in the current state of 
the system. is the sum of the rates of all nodes in the left subtree of node k and r^ is 
the sum of all rates in the right subtree of node k. Using these quantities, the sum of all 
rates is R{s) = wi + li + ri. Formally, we define Wk = Ik = = for all k > M. Then, 
h = ''JJ2k + hk + i^2k and = W2k+i + hk+i + f2k+i- The following algorithm selects a random 
event by climbing down the tree. The index k points to the current node. 

Algorithm B.2 Selection of events with probability Pk{s) = Wk/{wi + h + ri) 

1. Draw a uniformly distributed random number ^ G [0,wi + h + ri[. Set A; <— 1 (Start at 
the root node). 

2. If i < Ik, then set k <— 2k, i.e. go to the left child. Else, if < Ik + Wk return event 
k. In this case, the algorithm is done. Else, go to the right child and adjust ^ such that 
^<i<rk: k^2k + l, i-h-Wk. 
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3. Go to step 2. 



The proof that algorithm is correct is by induction. Whenever the algorithm enters node 
k in step 2, the actual value of ^ is a uniformly distributed random number G [0, u)fc + 4 + '"fcf- 
This is ensured by the initialization in step one and the adjustment in step 2. Therefore, the 
algorithm terminates with probability Wk/{wk + ^fc + fk)- The probabilities to go to the left 
and the right child are lk/{wk + Ik + fk) and rk/{wk + h + ^k), respectively. In particular, 
the probability to terminate at the root is wi/{wi + li + ri). This is the basis of induction. 
Node /c is a child of node m = [k/2\ < k. By induction, the algorithm terminates in node m 
with the correct probability, which can only be the case if the algorithm enters node m with 
probability (wm + Im + Tm)/{wi + li + ri). The probability to go from node m to node k is 
{wk + lk + fk)/{wm + Im + rm)- Therefore, the probability to terminate in node k is 

Wra + lm + Tm Wk + k + Tk Wk Wk , 

1 • ; • ; = ; • q-e.d. 

Wi+li+ n Wm + lm + rm Wk + k + Tk Wi + h + Ti 

In each iteration where the algorithm does not terminate, one half of the events in the subtrees 
of node k are ruled out. Therefore, the average running time of the algorithm is 0(log2 AA). 



In step 2 of algorithm B.l, the tree must be updated to reflect the changes in the Pi{s). 
To do this for one event, we set Wi = Pi{s). Then, the Ik and of the parents of node i must 
be adjusted. This is done by the following algorithm which, on average, requires C'(log2A^) 
iterations. 

Algorithm B.3 Update the rate of event i 

1. Start at node i. Set k ^ i, Wi ^ Pi{s)- 

2. If k = 1, terminate. Else, go to the parent node by setting k <— [k/2\ . 

3. Adjust Ik and rk- Set Ik ^ W2k + hk + r2k and rk ^ W2k+i + hk+i + ?'2fc+i- 
4- Go to step 2. 

If the rates of events change in each iteration, the time required for the complete update 
isO(log2AA). 

These algorithms work in the general case where there is a unique rate for every event. 
However, in many lattice gas models the rate of an event depends only on the local environment 
of the sites which are modified by the event. Then, the number of different rates M is 
independent of the system size and much smaller than the number of events N . This feature 
can be exploited to make the running time both of the random selection of events and the 
update of data structures 0{1). The idea is to keep lists of events the rates of which are equal. 
Consider the case where we have M lists {lg}^i. If the ng{s) events in list Ig have the rate 
Pg, the total rate of these events is Rg{s) = ng{s)pg. The probability for one member of the 
list to be selected is 

Pg ^ Rgjs) 1 
R{s) R{s) ' ng{s)- 

Consequently, we can perform the selection of events in a two-step process. First, we select a 
list Ig with probability Rg{s)/R{s). This can be done with a binary tree where each node is 



attributed to one list using algorithm B.2. Then, we select one of the members of L randomly 
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with equal probability for all events. On update, events the rates of which have changed 
are moved to their new group. Then, the Rg{s) of lists the length of which has changed are 



updated using algorithm B.3. Since the members of the lists can be accessed in 0(1) time, the 
running time both for the selection of events and the updates is independent of M. However, 
compared to the direct storage of the rates of events in a binary tree, the manipulation of the 
lists requires some computational overhead. In our experience, in a two-dimensional lattice 
gas at system sizes of a few hundred lattice constants, this overhead outweighs the gain of 
speed whenever Ai is greater than a few dozen. 

B.2 Application to lattice gases in thermal equilibrium 

The application to grand-canonical simulations of the lattice gases discussed in chapter ^ is 
straightforward. In a model where each site is in one of n possible states S {0, n — 1}, there 
are n — 1 events m G {1, ...,n — 1} per lattice site which change the state of the site: Xij <— 
(xjj -|- m) modn. This can be done with Metropolis-type rates Vg^s' = min{l, exp[—{H{s') — 
H(s))/{kT)]} or symmetrical rates r^^s' = exp[—(H(s') — H(s))/{2kT)]. Both possibilities 



fulfil the detailed balance condition [B.3| and yield identical results. 

A canonical simulation requires an algorithm where the number of particles is fixed. In 
this case, in every event one particle jumps to a different site. We choose a nonlocal dynamics 
where the range of particle jumps is unlimited. This yields considerably faster equilibration 
compared to a Kawasaki dynamics with nearest neighbour diffusion only. For simplicity, we 
permit only jumps to a site where the binding energy of the particle is independent of the 
state of its initial site, i.e. we forbid jumps to nearest and next nearest neighbour sites. If the 
particle jumps from site i to site j, the energy difference between the final and the initial state 
is AH = AHj — AHi, where AHx is the energy difference of the system with site x occupied 
and empty. The rates 

n^j = exp [{AHi - AHj) I {2kT)] (B.IO) 

fulfil the detailed balance condition. Then, the probability for a jump from site i to site j 
factorizes, i.e. 

^± 

Pi^3 = Pi ■ Pj '"'here = ^ ^ . 

Here, we have introduced the rate for deposition of a particle at site x (r^) and for removal 
of a particle at site x (r~). = exp[— Affa;/(2/cr)] if site x is empty and zero otherwise. 
Conversely, r~ = eyi.-p[AHx/ {"ikT)] on occupied sites and zero on empty sites. Due to this 
factorization property we can proceed in two steps: In the first step, we select the site i from 
which the particle starts with probability p^. Then, we select the site j where the particle 
is landing with probability p^. If the distance between site i and site j is > 2 the particle 
is moved. Otherwise, the event is rejected and the system remains unchanged. Since the 
number of rejected events is small on large systems, the loss of speed can be neglected. 

In our model of flat CdTe(OOl) surfaces only the number of Cd atoms is fixed, while the 
number of Te dimers may change. This requires a more elaborate algorithm which uses both 
canonical and grand-canonical techniques. A Cd atom may jump to any site which is not 
occupied by a Cd atom. If the arrival site is occupied by a dimer, the dimer is destroyed. 
The removal of the Cd atom at the starting site creates a pair of Te atoms. We consider 
both the case where these Te atoms dimerize immediately and the case where they remain 
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unbound. Additionally, Te dimers may break up and dimers may be created on empty sites. 
In the following, these processes will be denoted as "dimer flips". In each timestep, a random 
number is drawn to decide whether a Cd jump with immediate formation of a dimer at the 
starting site, a Cd jump without dimerization or a dimer flip will occur. The probabilities are 
proportional to the sums of the rates of the processes in each group. A Lx N system contains 
0{LN) Cd atoms and Te dimers. Thus, there are 0(LN) dimer flip events. Since there are 
0{LN) possible arrival sites for each Cd atom, there are 0{L?'N'^) Cd jump events. To keep 
the ratio of dimer flips and Cd jumps independent of the system size, the dimer flips have 
been weighted with a prefactor LN . Then, one event in the selected group is perfomed. For 
dimer flips, this is done with the grand-canonical algorithm while Cd jumps are performed by 
the canonical two-step algorithm. 



B.3 Application to kinetical Monte Carlo simulations 

In kinetical Monte Carlo simulations of processes like crystal growth, the set of events i G 
{1, A^} which can occur in a particular state of the system and their rates pi are determined 
by the physics of the investigated model. Typically, this set of events includes the deposition 
of adatoms at the surface, diffusion of particles at the surface and desorption. In all full 
diffusion Monte Carlo simulations which are presented in this thesis, we assume that the rates 



of diffusion and desorption are given by the Arrhenius law 1.2. For a physically reasonable 



choice of the activation energies, the rates for diffusion processes fulfil the detailed balance 



condition B.3 such that a system where there is neither adsorption nor desorption will finally 



reach thermal equilibrium. 
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